2009-01-04 25 views
5

这是继续发布的问题: Finding the center of mass on a 2D bitmap其中讨论了如何在布尔矩阵中找到质心,如给出的示例。查找矩阵/位图中的质量簇

假设现在我们的矩阵扩展到这种形式:

0 1 2 3 4 5 6 7 8 9 
1 . X X . . . . . . 
2 . X X X . . X . . 
3 . . . . . X X X . 
4 . . . . . . X . . 
5 . X X . . . . . . 
6 . X . . . . . . . 
7 . X . . . . . . . 
8 . . . . X X . . . 
9 . . . . X X . . . 

正如你可以看到我们现在有4个中心的质量,为4个不同的集群。

我们已经知道如何找到一个质量中心,因为只有一个存在,如果我们在这个矩阵上运行该算法,我们会在矩阵的中间得到一些点,这对我们没有帮助。

什么可以是一个很好,正确和快速的算法来找到这些质量群?

回答

3

我想我会检查矩阵中的每个点,并根据它的邻居找出它的质量。点的质量会随着距离的平方而下降。然后你可以选择距离彼此最小的前四个点。

下面是一些Python代码,我一起试着说明找出每个点的质量的方法。

matrix = [[1.0 if x == "X" else 0.0 for x in y] for y in """.XX...... 
.XXX..X.. 
.....XXX. 
......X.. 
.XX...... 
.X....... 
.X....... 
....XX... 
....XX...""".split("\n")] 

HEIGHT = len(matrix) 
WIDTH = len(matrix[0]) 
Y_RADIUS = HEIGHT/2 
X_RADIUS = WIDTH/2 

要计算质量对于给定的点:

def distance(x1, y1, x2, y2): 
    'Manhattan distance http://en.wikipedia.org/wiki/Manhattan_distance' 
    return abs(y1 - y2) + abs(x1 - x2) 

def mass(m, x, y): 
    _mass = m[y][x] 
    for _y in range(max(0, y - Y_RADIUS), min(HEIGHT, y + Y_RADIUS)): 
    for _x in range(max(0, x - X_RADIUS), min(WIDTH, x + X_RADIUS)): 
     d = max(1, distance(x, y, _x, _y)) 
     _mass += m[_y][_x]/(d * d) 
    return _mass 

注:我使用Manhattan距离(又名Cityblock,又名出租车几何)在这里,因为我不使用你的榜样矩阵一些设置认为使用欧几里德距离的附加精度值得调用sqrt()的代价。

迭代通过我们的矩阵和建立一个元组列表像(X,Y,质量(X,Y)):

point_mass = [] 
for y in range(0, HEIGHT): 
    for x in range(0, WIDTH): 
    point_mass.append((x, y, mass(matrix, x, y))) 

排序的质量对每个点的地址列表:

from operator import itemgetter 
point_mass.sort(key=itemgetter(2), reverse=True) 

望着在排序列表顶部9点:

(6, 2, 6.1580555555555554) 
(2, 1, 5.4861111111111107) 
(1, 1, 4.6736111111111107) 
(1, 4, 4.5938888888888885) 
(2, 0, 4.54) 
(4, 7, 4.4480555555555554) 
(1, 5, 4.4480555555555554) 
(5, 7, 4.4059637188208614) 
(4, 8, 4.3659637188208613) 

如果我们将努力从最高到最低,并过滤接[R远点过于接近已经看到了,我们会得到点(我做手工,因为我已经用完了时间,现在做在代码中...):

(6, 2, 6.1580555555555554) 
(2, 1, 5.4861111111111107) 
(1, 4, 4.5938888888888885) 
(4, 7, 4.4480555555555554) 

哪个一个非常直观的结果,只需查看矩阵(请注意,与您的示例进行比较时,坐标为零)。

1

我首先想到的是首先找到任何非零值的单元格。从那里做一些洪水填充算法,并计算找到的细胞的质心。接下来从矩阵中找出找到的单元格,并从顶部开始。

这当然不会像Google的方法那样扩展,即tuinstoel链接,但对于较小的矩阵更容易实现。

编辑:

Disjoint sets(使用路径压缩和工会通过秩)也许会有帮助。它们具有øαÑ))时为联合复杂性和查找集,其中

αÑ)= MIN {ķ:甲ķ(1)≥ n}。

ķÑ)是阿克曼功能,所以αÑ)将基本上是ø(1)对于任何合理的值。唯一的问题是不相交集合是项目到集合的单向映射,但是如果你要通过所有项目,这并不重要。

这里是一个演示简单的Python脚本:

from collections import defaultdict 

class DisjointSets(object): 
    def __init__(self): 
     self.item_map = defaultdict(DisjointNode) 

    def add(self,item): 
     """Add item to the forest.""" 
     # It's gets initialized to a new node when 
     # trying to access a non-existant item. 
     return self.item_map[item] 

    def __contains__(self,item): 
     return (item in self.item_map) 

    def __getitem__(self,item): 
     if item not in self: 
      raise KeyError 
     return self.item_map[item] 

    def __delitem__(self,item): 
     del self.item_map[item] 

    def __iter__(self): 
     # sort all items into real sets 
     all_sets = defaultdict(set) 
     for item,node in self.item_map.iteritems(): 
      all_sets[node.find_set()].add(item) 
     return all_sets.itervalues() 

class DisjointNode(object): 
    def __init__(self,parent=None,rank=0): 
     if parent is None: 
      self.parent = self 
     else: 
      self.parent = parent 
     self.rank = rank 

    def union(self,other): 
     """Join two sets.""" 
     node1 = self.find_set() 
     node2 = other.find_set() 
     # union by rank 
     if node1.rank > node2.rank: 
      node2.parent = node1 
     else: 
      node1.parent = node2 
      if node1.rank == node2.rank: 
       node2.rank += 1 
     return node1 

    def find_set(self): 
     """Finds the root node of this set.""" 
     node = self 
     while node is not node.parent: 
      node = node.parent 
     # path compression 
     root, node = node, self 
     while node is not node.parent: 
      node, node.parent = node.parent, root 
     return root 

def find_clusters(grid): 
    disj = DisjointSets() 
    for y,row in enumerate(grid): 
     for x,cell in enumerate(row): 
      if cell: 
       node = disj.add((x,y)) 
       for dx,dy in ((-1,0),(-1,-1),(0,-1),(1,-1)): 
        if (x+dx,y+dy) in disj: 
         node.union(disj[x+dx,y+dy]) 
    for index,set_ in enumerate(disj): 
     sum_x, sum_y, count = 0, 0, 0 
     for x,y in set_: 
      sum_x += x 
      sum_y += y 
      count += 1 
     yield 1.0 * sum_x/count, 1.0 * sum_y/count 

def main(): 
    grid = [[('.' != cell) for cell in row if not cell.isspace()] for row in (
     ". X X . . . . . .", 
     ". X X X . . X . .", 
     ". . . . . X X X .", 
     ". . . . . . X . .", 
     ". X X . . . . . .", 
     ". X . . . . . . .", 
     ". X . . . . . . .", 
     ". . . . X X . . .", 
     ". . . . X X . . .", 
    )] 
    coordinates = list(find_clusters(grid)) 
    centers = dict(((round(x),round(y)),i) for i,(x,y) in enumerate(coordinates)) 
    for y,row in enumerate(grid): 
     for x,cell in enumerate(row): 
      if (x,y) in centers: 
       print centers[x,y]+1, 
      elif cell: 
       print 'X', 
      else: 
       print '.', 
     print 
    print 
    print '%4s | %7s %7s' % ('i','x','y') 
    print '-'*22 
    for i,(x,y) in enumerate(coordinates): 
     print '%4d | %7.4f %7.4f' % (i+1,x,y) 

if __name__ == '__main__': 
    main() 

输出:

. X X . . . . . . 
. X 3 X . . X . . 
. . . . . X 4 X . 
. . . . . . X . . 
. X X . . . . . . 
. 2 . . . . . . . 
. X . . . . . . . 
. . . . X X . . . 
. . . . X 1 . . . 

    i |  x  y 
---------------------- 
    1 | 4.5000 7.5000 
    2 | 1.2500 4.7500 
    3 | 1.8000 0.6000 
    4 | 6.0000 2.0000 

的这点是要证明不相交的集合。 find_clusters()中的实际算法可以升级为更强大的功能。

参考

  • 算法导论。 2nd ed。 Cormen et.al.
1

Here's一个类似的问题,一个不那么快的算法,以及其他几个更好的方法来做到这一点。