Procrustes Análisis con NumPy?

¿Hay algo como la función de procrustes de Matlab en NumPy / SciPy o bibliotecas relacionadas?


Para referencia. El análisis de Procrustes apunta a alinear 2 conjuntos de puntos (en otras palabras, 2 formas) para minimizar la distancia cuadrada entre ellos mediante la eliminación de los componentes de distorsión de escala, traslación y rotación.

Ejemplo en Matlab:

 X = [0 1; 2 3; 4 5; 6 7; 8 9]; % first shape R = [1 2; 2 1]; % rotation matrix t = [3 5]; % translation vector Y = X * R + repmat(t, 5, 1); % warped shape, no scale and no distortion [d Z] = procrustes(X, Y); % Z is Y aligned back to X Z Z = 0.0000 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 

La misma tarea en NumPy:

 X = arange(10).reshape((5, 2)) R = array([[1, 2], [2, 1]]) t = array([3, 5]) Y = dot(X, R) + t Z = ??? 

Nota: Solo me interesa la forma alineada, ya que el error cuadrado (variable d en el código de Matlab) se calcula fácilmente a partir de 2 formas.

No tengo conocimiento de ninguna implementación preexistente en Python, pero es fácil echar un vistazo al código MATLAB usando edit procrustes.m y trasladarlo a Numpy:

 def procrustes(X, Y, scaling=True, reflection='best'): """ A port of MATLAB's `procrustes` function to Numpy. Procrustes analysis determines a linear transformation (translation, reflection, orthogonal rotation and scaling) of the points in Y to best conform them to the points in matrix X, using the sum of squared errors as the goodness of fit criterion. d, Z, [tform] = procrustes(X, Y) Inputs: ------------ X, Y matrices of target and input coordinates. they must have equal numbers of points (rows), but Y may have fewer dimensions (columns) than X. scaling if False, the scaling component of the transformation is forced to 1 reflection if 'best' (default), the transformation solution may or may not include a reflection component, depending on which fits the data best. setting reflection to True or False forces a solution with reflection or no reflection respectively. Outputs ------------ d the residual sum of squared errors, normalized according to a measure of the scale of X, ((X - X.mean(0))**2).sum() Z the matrix of transformed Y-values tform a dict specifying the rotation, translation and scaling that maps X --> Y """ n,m = X.shape ny,my = Y.shape muX = X.mean(0) muY = Y.mean(0) X0 = X - muX Y0 = Y - muY ssX = (X0**2.).sum() ssY = (Y0**2.).sum() # centred Frobenius norm normX = np.sqrt(ssX) normY = np.sqrt(ssY) # scale to equal (unit) norm X0 /= normX Y0 /= normY if my < m: Y0 = np.concatenate((Y0, np.zeros(n, m-my)),0) # optimum rotation matrix of Y A = np.dot(X0.T, Y0) U,s,Vt = np.linalg.svd(A,full_matrices=False) V = Vt.T T = np.dot(V, UT) if reflection is not 'best': # does the current solution use a reflection? have_reflection = np.linalg.det(T) < 0 # if that's not what was specified, force another reflection if reflection != have_reflection: V[:,-1] *= -1 s[-1] *= -1 T = np.dot(V, UT) traceTA = s.sum() if scaling: # optimum scaling of Y b = traceTA * normX / normY # standarised distance between X and b*Y*T + c d = 1 - traceTA**2 # transformed coords Z = normX*traceTA*np.dot(Y0, T) + muX else: b = 1 d = 1 + ssY/ssX - 2 * traceTA * normY / normX Z = normY*np.dot(Y0, T) + muX # transformation matrix if my < m: T = T[:my,:] c = muX - b*np.dot(muY, T) #transformation values tform = {'rotation':T, 'scale':b, 'translation':c} return d, Z, tform 

Hay una función de Scipy para ello: scipy.spatial.procrustes

Solo estoy publicando su ejemplo aquí:

 >>> import numpy as np >>> from scipy.spatial import procrustes >>> a = np.array([[1, 3], [1, 2], [1, 1], [2, 1]], 'd') >>> b = np.array([[4, -2], [4, -4], [4, -6], [2, -6]], 'd') >>> mtx1, mtx2, disparity = procrustes(a, b) >>> round(disparity) 0.0