结构化3D数据聚类

Question

假设我有一些物体（类似于蛋白质，但不完全相同），每个物体都用一组 n 个三维坐标来表示。这些物体在空间中有不同的方向。我们可以通过使用Kabsch算法来对齐它们，然后计算对齐后坐标的均方根偏差，从而判断它们的相似性。

我想知道，如何对一大堆这样的结构进行聚类，以找出最常见的那一类（也就是大多数结构所属的那一类）。另外，有没有办法在 Python 中实现这个？举个简单的例子，这里有一组没有聚类的结构（每个结构由四个顶点的坐标表示）：

然后是我想要的聚类结果（使用两个聚类）：

我试过把所有结构对齐到一个参考结构（也就是第一个结构），然后使用 Pycluster.kcluster 在参考结构和对齐后坐标之间的距离上进行 k-means 聚类，但这似乎有点笨拙，效果也不太好。每个聚类中的结构彼此之间并不太相似。理想情况下，这种聚类应该是在实际结构上进行，而不是在差异向量上，但这些结构的维度是 (n,3)，而 k-means 聚类需要的是 (n,) 的维度。

我尝试的另一个选项是 scipy.clustering.hierarchical。这个方法似乎效果不错，但我在决定哪个聚类是最常见的时遇到了麻烦，因为总是可以通过向上移动到树的下一个分支找到一个更常见的聚类。

如果有任何想法、建议或关于不同（已经在 Python 中实现的）聚类算法的点子，我将非常感激。

Answer 1

为了给我自己提问的内容一个简单的介绍，我想说可以用形状中每个点之间的距离列表作为聚类的标准。

我们先来创建一些形状：

shapes = np.array([[[1,4],[4,2],[11,2],[14,0]],
          [[4,5],[4,2],[11,2],[13,0]],
          [[1,3],[4,2],[11,2],[14,1.5]],
          [[3,5],[4,2],[10,7],[7,9]],
          [[5,5],[4,2],[10,7],[6,6]]])

def random_rotation():
    theta = 3 * np.pi * np.random.random()
    rotMatrix = numpy.array([[np.cos(theta), -np.sin(theta)], 
                             [np.sin(theta),  np.cos(theta)]])
    return rotMatrix

new_shapes = []
for s in shapes:
    rr = random_rotation()
    new_shapes += [[list(rr.dot(p)) + [0] for p in s]]
new_shapes = np.array(new_shapes)

for i, s in enumerate(new_shapes):
    plot(s[:,0], s[:,1], color='black')
    text(np.mean(s[:,0]), np.mean(s[:,1]), str(i), fontsize=14)

接着，我们创建一些辅助函数，并生成一个数组，里面包含每个形状中所有顶点之间的距离（darray）。

import itertools as it

def vec_distance(v1, v2):
    '''
    The distance between two vectors.
    '''
    diff = v2 - v1
    return math.sqrt(sum(diff * diff))

def distances(s):
    '''
    Compute the distance array for a shape s.
    '''
    ds = [vec_distance(p1, p2) for p1,p2 in it.combinations(s, r=2)]

    return np.array(ds)


# create an array of inter-shape distances for each shape
darray = np.array([distances(s) for s in new_shapes])

然后使用 Pycluster 将它们分成两个聚类。

import Pycluster as pc

clust = pc.kcluster(darray,2)
print clust

结果发现，第一个聚类里有三个形状，第二个聚类里有两个形状。

(array([0, 0, 0, 1, 1], dtype=int32), 4.576996142441375, 1)

那么这些形状分别是什么呢？

import brewer2mpl

dark2 = brewer2mpl.get_map('Dark2', 'qualitative', 4).mpl_colors

for i,(s,c) in enumerate(zip(new_shapes, clust[0])):
    plot(s[:,0], s[:,1], color=dark2[c])
    text(np.mean(s[:,0]), np.mean(s[:,1]), str(i), fontsize=14)

看起来不错！不过问题是，随着形状变大，距离数组的大小会以顶点数量的平方速度增长。我找到了一份演示文稿，里面描述了这个问题，并提出了一些解决方案（比如我猜是某种降维方法的SVD）来加快处理速度。

我还不打算接受我的答案，因为我想听听其他人对这个简单问题的想法或建议。