Grassmann module

Grassmann.Karcher(shapes, eps=1e-08, max_steps=20)

Karcher mean for given shapes.

Calculated Karcher mean for given shapes (elements on Grassmann) by minimizing the sum of squared (Riemannian) distances to all shapes in the data (Fletcher, Lu, and Joshi 2003)

Parameters:

  • shapes – (n_shapes, n_landmarks, 2) array of given shapes (Grassmann elements)

  • max_steps – maximum number of iterations to converge

Returns:

(n_landmarks, 2) array defining Karcher mean (element on Grassmann)

Grassmann.PGA(mu, shapes_gr, n_coord=None)

Principal Geodesic Analysis (PGA).

Principal Geodesic Analysis (PGA), a generalization of Principal Component Analysis (PCA) over Riemannian manifolds. PGA is a data-driven approach that determines principal components as elements in a central tangent space, given a data set represented as elements in asmooth manifold.

Parameters:

  • mu – (n_landmarks, 2) array defining Karcher mean (element on Grassmann)

  • shapes_gr – (n_shapes, n_landmarks, 2) given shapes (elements on Grassmann)

  • n_coord – dimension of resulting PGA space (if None n_coord=n_landmarks)

Returns:

Vh is principal basis transposed ((n_coord*2 - 4)x(n_coord*2 - 4)), t are given elements in principal coordinates, S is corresponding singular values,

Grassmann.PGA_modes(PGA_directions, mu, scale=0.1, sub=10)

Moves given element on Grassmann in each of given directions.

Parameters:

  • PGA_directions – directions to perturb element mu

  • scale – how much we perturb element mu

  • mu – element on grassmann according which we perturb

  • sub – subset of directions (first sub directions)

Returns:

elements on Grassmann (perturbed from mu in PGA_directions)

Grassmann.check_orthogonality(X)

Checks orthogonality of a matrix.

Parameters:

X – np.array of matrix

Raises:

ValueError – raises error is given matrix is not orthogonal

Grassmann.distance(X, Y)

Geodesic distance on Grassmannian.

Geodesic distance on Grassmannian is defined by the principal angles between the subspaces spanned by the columns of X and Y, denoted by span(X) and span(Y). The cosines of the principal angles theta1 and theta2 between span(X) and snap(Y) are the singular values of X.T@Y. That is, X.T@Y = U D V.T, where D = diag(cos(theta1), cos(theta2)). The distance between two shapes is then defined as dist = sqrt(theta1**2+theta2**2).

Parameters:

  • X – (n_landmarks, 2) array defining first shape

  • Y – (n_landmarks, 2) array defining second shape

Returns:

distance between two shapes on Grassmannian

Grassmann.exp(t, X, direction)

Exponential mapping (Grassmannian geodesic).

Parameters:

  • t – scalar > 0, how far in given direction to move (if t=0, exp(X, log_map) = X)

  • X – (n_landmarks, 2) array defining starting point of geodesic on Grassmann

  • direction – (n_landmarks, 2) array defining direction in tangent space (tangent vector Delta)

Returns:

(n_landmarks, 2) array defining end point on Grassmann

Grassmann.get_PGA_coordinates(shapes_gr, mu, V, n_modes=None)

Get PGA coordinates of given standardized shapes (elements of Grassmann) and PGA space defined by Karcher mean and basis vectors.

Parameters:

  • shapes_gr – (n_shapes, n_landmarks, 2) array of given standardized shapes (element of Grassmann)

  • mu – (n_landmarks, 2) array of Karcher mean (origin of PGA space)

  • V – (n_coord*2, n_coord*2) array of PGA basis vectors

Returns:

(n_shapes, n_coord) array of PGA coordinates for given shapes

Grassmann.landmark_affine_transform(X_phys)

Shift and scale all shapes using Landmark-Affine standardization (Bryner, 2D affine and projective spaces).

LA-standardization normalizes the shape such that it has zero mean (translation invariance) and sample covariance proportional to I2 over the n discrete boundary landmarks defining the shape.

Parameters:

X_phys – (n_shapes, n_landmarks, 2) array of physical coordinates defining shapes

Returns:

X_grassmann, M, b, such that X_phys = X_grassmann @ M + b.

Grassmann.log(X, Y)

Logarithmic mapping (inverse mapping of exponential map). Algorithm 11 (Zimmermann, 2019)

Calculate logarithmic map log_X(Y) (inverse mapping of exponential map). Calculates direction(tangent vector Delta) from X to Y in tangent subspace.

Parameters:

  • X – (n_landmarks, 2) array defining start point of geodesic on Grassmann

  • Y – (n_landmarks, 2) array defining end point of geodesic on Grassmann

Returns:

(n_landmarks, 2) array defining direction in tangent space (tangent vector Delta)

Grassmann.parallel_translate(start, end_direction, vector)

Parallel translation of a vector along geodesic from start point to end point Edelman et al. (Theorem 2.4, pp. 321).

Parameters:

  • start – (n_landmarks, 2) array defining start point element on Grassmann

  • end_direction – (n_landmarks, 2) array defining direction to the end point element on Grassmann (log map)

  • vector – (n_landmarks, 2) array defining vector to be translated

Returns:

(n_landmarks, 2) array defining translated vector

Grassmann.perturb_gr_shape(Vh, mu, perturbation)

Given element Karcher mean, perturbs it in given direction by a given amount.

Parameters:

  • Vh – (n_landmarks*2 - 4 , n_landmarks*2 - 4) array of PGA basis vectors transposed

  • mu – (n_landmarks, 2) array of Karcher mean (elenemt on Grassmann)

  • perturbation – (n_modes,) array of amount of perturbations in pga coordinates

Returns:

(n_landmarks, 2) array of perturbed element on Grassmann

Grassmann.procrustes(X, Y)

Procrustes clustering match the shapes via Procrustes analysis (Gower 1975).

This function calculates rotation that can be applied to shapes for matching (rotation does not fundamentally modify the elements in the Grassmannian).

Parameters:

  • X – (n_landmarks, 2) array defining shape 1

  • Y – (n_landmarks, 2) array defining shape 2

Returns:

(2, 2) array of rotation matrix