Анализ прокрустов с помощью NumPy?

Question

Есть ли что-то вроде функции Matlab procrustes в NumPy/SciPy или связанных с ней библиотеках?

---

Для справки. Анализ Procrustes направлен на выравнивание двух наборов точек (другими словами, 2 фигуры), чтобы минимизировать квадратное расстояние между ними, удаляя компоненты масштабирования, трансляции и вращения.

Пример в 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 

же задача в NumPy:

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

Примечание: Я заинтересован только в соосных форме, так как ошибки (переменная d в Matlab коде) легко вычисленный из двух форм.

Answer 1

Я не знаю ни уже существующей реализации в Python, но это легко взглянуть на код MATLAB, используя edit procrustes.m и портирование на 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, U.T) 

    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, U.T) 

    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 

Answer 2

Существует функция SciPy для это: scipy.spatial.procrustes

Я просто разместить свой пример здесь:

>>> 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