多维标度法(MDS,Multidimensional Scaling)及普氏分析(Procrustes Analysis)在人体姿态关节点上的简单示例(python)

多维标度法(MDS,Multidimensional Scaling)

多维标度法一个简单的应用示例就是，已知一组城市之间的相对距离关系(相似矩阵)，如何求解出各个城市在地图上的位置，使其尽可能满足前面的相对距离约束。有关多维标度法的具体介绍大家可以参考网上其他博主的讲解，我在这里就不具体展开了。

普氏分析(Procrustes Analysis)

最近网上比较火的换脸操作大部分都会用到普氏分析法，它是一种用来分析形状分布的方法，寻找标准形状，利用最小二乘法寻找样本形状到标准形状的仿射变化。举换脸的例子就是，两张人脸上有各自的landmarks，现在要讲B脸和A脸对齐(对应的landmark要匹配)，这样换得的结果才更加真实。但是由于A B两张人脸的大小，旋转等不太相同，如何变换呢？答案是使用普氏分析法寻找B到A的一个变换。

在人体姿态关节点上的简单示例(python实现，图片demo在最下方)

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets, decomposition, manifold
# from scipy.spatial import procrustes

def main():
	# human3.6m的关节连接顺序
    human36m_connectivity_dict = [(0, 1), (1, 2), (2, 6), (5, 4), (4, 3), (3, 6), (6, 7), (7, 8), (8, 16), (9, 16),
                                  (8, 12),
                                  (11, 12), (10, 11), (8, 13), (13, 14), (14, 15)]
	
	# 定义figure以及子图
    fig = plt.figure(figsize=(10, 10))
    ax1 = fig.add_subplot(2, 2, 1)
    ax1.set_title('original') # 原始骨架
    ax2 = fig.add_subplot(2, 2, 2)
    ax2.set_title('MDS reconstructed') # MDS重构出的骨架
    ax3 = fig.add_subplot(2, 2, 3)
    ax3.set_title('procrustes') # Procrustes对齐之后的骨架
    ax4 = fig.add_subplot(2, 2, 4)
    ax4.set_title('alignment of original and procrustes') # MDS重构出的骨架和Procrustes变换之后的骨架进行对齐

	# 定义的一组2d关节点坐标（human3.6m， shape：17 x 2）
    keypoints = 1000 - np.array([[524, 491],
                                 [521, 437],
                                 [522, 378],
                                 [548, 385],
                                 [552, 446],
                                 [561, 502],
                                 [535, 382],
                                 [535, 348],
                                 [537, 315],
                                 [537, 286],
                                 [462, 323],
                                 [490, 321],
                                 [519, 315],
                                 [560, 318],
                                 [594, 328],
                                 [618, 334],
                                 [533, 300]])

    for (index_from, index_to) in human36m_connectivity_dict:
        xs, ys = [np.array([keypoints[index_from, j], keypoints[index_to, j]]) for j in range(2)]
        ax1.plot(xs, ys, linewidth=2.5)
	
    mds = manifold.MDS(n_components=2)
    # MDS求解 X_r即为求解得到的关节点的坐标，其尽可能的满足了原始骨架坐标之间的距离约束，mds.fit_transform(keypoints)内部会根据keypoints自动生成一个距离相似矩阵，默认为欧式距离，后面的fit过程便是在这个相似距离矩阵上进行的。
    X_r = mds.fit_transform(keypoints)

    for (index_from, index_to) in human36m_connectivity_dict:
        xs, ys = [np.array([X_r[index_from, j], X_r[index_to, j]]) for j in range(2)]
        ax2.plot(xs, ys, linewidth=2.5)
	
	# Procrustes 将X_r与原始keypoints对齐，Z_pts为X_r对齐之后的坐标。
    d, Z_pts, Tform = procrustes(keypoints, X_r)

    for (index_from, index_to) in human36m_connectivity_dict:
        xs, ys = [np.array([Z_pts[index_from, j], Z_pts[index_to, j]]) for j in range(2)]
        ax3.plot(xs, ys, linewidth=2.5)

    for (index_from, index_to) in human36m_connectivity_dict:
        xs, ys = [np.array([Z_pts[index_from, j], Z_pts[index_to, j]]) for j in range(2)]
        ax4.plot(xs, ys, linewidth=2.5)
    # keypoints[:, 0] = keypoints[:, 0] + 10
    for (index_from, index_to) in human36m_connectivity_dict:
        xs, ys = [np.array([keypoints[index_from, j], keypoints[index_to, j]]) for j in range(2)]
        ax4.plot(xs, ys, linewidth=2.5)

    plt.show()


def procrustes(X, Y, scaling=True, reflection='best'):
    """
       该函数是matlab版本对应的numpy实现
       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


if __name__ == '__main__':
    main()

由下图大家可以看到，MDS 重构出来的姿态和原始姿态大体相同，但是尺度和朝向会有差异，每次运行重构结果都不太一样。经过Procrustes Analysis对齐之后，尺度和朝向基本一致，但是也会有少许误差(该例子中很小，上述代码中procrustes函数返回的d值即为MSE误差)。