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

Есть ли что-то вроде функции 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) легко вычисляется из двух форм.

Ответ 1

Мне не известно о какой-либо ранее существовавшей реализации в Python, но легко взглянуть на код MATLAB с помощью edit procrustes.m и перенести его на Numpy:

def procrustes(X, Y, scaling=True, reflection='best'):
    """
    A port of MATLAB `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 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

Ответ 2

Для этого есть функция scipy.spatial.procrustes: 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

Ответ 3

scipy.spatial.prosrustes дает сходство. Если вы заинтересованы в матрице преобразования, взгляните на scipy.linalg.orthogonal_procrustes

https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.orthogonal_procrustes.html

> from scipy.linalg import orthogonal_procrustes
> A = np.array([[ 2,  0,  1], [-2,  0,  0]])
> R, sca = orthogonal_procrustes(A, np.fliplr(A))
> R
array([[-5.34384992e-17,  0.00000000e+00,  1.00000000e+00],
   [ 0.00000000e+00,  1.00000000e+00,  0.00000000e+00],
   [ 1.00000000e+00,  0.00000000e+00, -7.85941422e-17]])
> sca
9.0