【问题标题】：Iterate over 2d array in an expanding circular spiral在扩展的圆形螺旋中迭代二维数组
【发布时间】：2012-02-17 05:59:59
【问题描述】：

给定一个n by n 矩阵M，在行i 和列j，我想以圆形螺旋遍历所有相邻值。

这样做的目的是测试某个函数f，它依赖于M，以找到远离(i, j) 的半径，其中f 返回True。所以，f 看起来像这样：

def f(x, y):
    """do stuff with x and y, and return a bool"""

并且会这样调用：

R = numpy.zeros(M.shape, dtype=numpy.int)
# for (i, j) in M
for (radius, (cx, cy)) in circle_around(i, j):
    if not f(M[i][j], M[cx][cy]):
       R[cx][cy] = radius - 1
       break

其中circle_around 是返回（迭代器到）圆形螺旋中的索引的函数。因此，对于M 中的每个点，此代码将计算并存储从f 返回True 的那个点的半径。

如果有更有效的计算 R 的方法，我也愿意接受。

更新：

感谢所有提交答案的人。我编写了一个简短的函数来绘制 circle_around 迭代器的输出，以显示它们的作用。如果您更新答案或发布新答案，则可以使用此代码来验证您的解决方案。

from matplotlib import pyplot as plt
def plot(g, name):
    plt.axis([-10, 10, -10, 10])
    ax = plt.gca()
    ax.yaxis.grid(color='gray')
    ax.xaxis.grid(color='gray')

    X, Y = [], []
    for i in xrange(100):
        (r, (x, y)) = g.next()
        X.append(x)
        Y.append(y)
        print "%d: radius %d" % (i, r)

    plt.plot(X, Y, 'r-', linewidth=2.0)
    plt.title(name)
    plt.savefig(name + ".png")

结果如下： plot(circle_around(0, 0), "F.J"):

plot(circle_around(0, 0, 10), "WolframH"):

我将 Magnesium 的建议编码如下：

def circle_around_magnesium(x, y):
    import math
    theta = 0
    dtheta = math.pi / 32.0
    a, b = (0, 1) # are there better params to use here?
    spiral = lambda theta : a + b*theta
    lastX, lastY = (x, y)
    while True:
        r = spiral(theta)
        X = r * math.cos(theta)
        Y = r * math.sin(theta)
        if round(X) != lastX or round(Y) != lastY:
            lastX, lastY = round(X), round(Y)
            yield (r, (lastX, lastY))
        theta += dtheta

plot(circle_around(0, 0, 10), "magnesium"):

如您所见，满足我正在寻找的接口的结果都没有产生一个圆形螺旋，该螺旋覆盖了 0、0 附近的所有索引。FJ 是最接近的，尽管 WolframH 的命中点是正确的，只是不是按螺旋顺序。

【问题讨论】：

您能否确认您的数组非常大，或者您必须多次这样做，或者真值函数很昂贵或者......？我可以想出一些东西，但这似乎是过早的优化，除非你真的需要避免在半径之外进行测试。当然，简单的解决方案是找到数组中每个假点的半径，然后找到最小的半径。如果你真的需要它的话，这是一个很好的问题。
@yakiimo，数组有 1-2 百万个条目。
F.J 的答案对你有用还是你需要一个真正的圈子？
我更喜欢真正的圆圈。
仅供参考，我的待办事项清单上还有这个，但现在很长。

标签： python matrix loops geometry spiral

【解决方案1】：

由于提到点的顺序无关紧要，我只是按照它们在给定半径处出现的角度 (arctan2) 对它们进行排序。更改N 以获得更多积分。

from numpy import *
N = 8

# Find the unique distances
X,Y = meshgrid(arange(N),arange(N))
G = sqrt(X**2+Y**2)
U = unique(G)

# Identify these coordinates
blocks = [[pair for pair in zip(*where(G==idx))] for idx in U if idx<N/2]

# Permute along the different orthogonal directions
directions = array([[1,1],[-1,1],[1,-1],[-1,-1]])

all_R = []
for b in blocks:
    R = set()
    for item in b:
        for x in item*directions:
            R.add(tuple(x))

    R = array(list(R))

    # Sort by angle
    T = array([arctan2(*x) for x in R])
    R = R[argsort(T)]
    all_R.append(R)

# Display the output
from pylab import *
colors = ['r','k','b','y','g']*10
for c,R in zip(colors,all_R):
    X,Y = map(list,zip(*R))

    # Connect last point
    X = X + [X[0],]
    Y = Y + [Y[0],]
    scatter(X,Y,c=c,s=150)
    plot(X,Y,color=c)

axis('equal')
show()

给N=8：

更多积分N=16（对不起色盲）：

这显然接近一个圆并按半径增加的顺序击中每个网格点。

【讨论】：

【解决方案2】：

随着距离的增加产生点的一种方法是将其分解为简单的部分，然后合并这些部分的结果.很明显itertools.merge 应该进行合并。 简单的部分是列，因为对于固定的 x，点 (x, y) 可以仅通过查看 y 的值来排序。

以下是该算法的（简单）实现。请注意，使用平方欧几里得距离，并且包括中心点。最重要的是，只考虑 range(x_end) 中带有 x 的点 (x, y)，但我认为这对于您的用例来说是可以的（在上面的符号中，x_end 将是 n）。

from heapq import merge
from itertools import count

def distance_column(x0, x, y0):
    dist_x = (x - x0) ** 2
    yield dist_x, (x, y0)
    for dy in count(1):
        dist = dist_x + dy ** 2
        yield dist, (x, y0 + dy)
        yield dist, (x, y0 - dy)

def circle_around(x0, y0, end_x):
    for dist_point in merge(*(distance_column(x0, x, y0) for x in range(end_x))):
        yield dist_point

编辑：测试代码：

def show(circle):
    d = dict((p, i) for i, (dist, p) in enumerate(circle))
    max_x = max(p[0] for p in d) + 1
    max_y = max(p[1] for p in d) + 1
    return "\n".join(" ".join("%3d" % d[x, y] if (x, y) in d else "   " for x in range(max_x + 1)) for y in range(max_y + 1))

import itertools
print(show(itertools.islice(circle_around(5, 5, 11), 101)))

测试结果（点按circle_around产生的顺序编号）：

             92  84  75  86  94                
     98  73  64  52  47  54  66  77 100        
     71  58  40  32  27  34  42  60  79        
 90  62  38  22  16  11  18  24  44  68  96    
 82  50  30  14   6   3   8  20  36  56  88    
 69  45  25   9   1   0   4  12  28  48  80    
 81  49  29  13   5   2   7  19  35  55  87    
 89  61  37  21  15  10  17  23  43  67  95    
     70  57  39  31  26  33  41  59  78        
     97  72  63  51  46  53  65  76  99        
             91  83  74  85  93

编辑2：如果您确实需要i 的负值，请将cirlce_around 函数中的range(end_x) 替换为range(-end_x, end_x)。

【讨论】：

这看起来不像是螺旋形——请参阅我对问题的更新。我错过了什么？产生的是这样的：i.stack.imgur.com/h0mNa.png
@JasonSundram：1）我假设非负指数，因为你写了一个 n x n 矩阵。 2）我还假设相同距离的点的顺序并不重要。这就是我从你的问题和cmets中得到的印象。有错吗？
明白了。感谢您的澄清和代码，以显示您的答案是如何工作的。我还在帖子中为您的代码重新生成了图像，以更准确地描述它的工作原理。

【解决方案3】：

如果您遵循 x 和 y 螺旋索引，您会注意到它们都可以以递归方式定义。因此，很容易想出一个递归生成正确索引的函数：

def helicalIndices(n):
    num = 0
    curr_x, dir_x, lim_x, curr_num_lim_x = 0, 1, 1, 2
    curr_y, dir_y, lim_y, curr_num_lim_y = -1, 1, 1, 3
    curr_rep_at_lim_x, up_x = 0, 1
    curr_rep_at_lim_y, up_y = 0, 1

    while num < n:
        if curr_x != lim_x:
            curr_x +=  dir_x
        else:
            curr_rep_at_lim_x += 1
            if curr_rep_at_lim_x == curr_num_lim_x - 1:
                if lim_x < 0:
                    lim_x = (-lim_x) + 1
                else:
                    lim_x = -lim_x
                curr_rep_at_lim_x = 0
                curr_num_lim_x += 1
                dir_x = -dir_x
        if curr_y != lim_y:
            curr_y = curr_y + dir_y
        else:
            curr_rep_at_lim_y += 1
            if curr_rep_at_lim_y == curr_num_lim_y - 1:
                if lim_y < 0:
                    lim_y = (-lim_y) + 1
                else:
                    lim_y = -lim_y
                curr_rep_at_lim_y = 0
                curr_num_lim_y += 1
                dir_y = -dir_y
        r = math.sqrt(curr_x*curr_x + curr_y*curr_y)        
        yield (r, (curr_x, curr_y))
        num += 1

    hi = helicalIndices(101)
    plot(hi, "helicalIndices")

从上图可以看出，这正是所要求的。

【讨论】：

【解决方案4】：

这是circle_around() 的基于循环的实现：

def circle_around(x, y):
    r = 1
    i, j = x-1, y-1
    while True:
        while i < x+r:
            i += 1
            yield r, (i, j)
        while j < y+r:
            j += 1
            yield r, (i, j)
        while i > x-r:
            i -= 1
            yield r, (i, j)
        while j > y-r:
            j -= 1
            yield r, (i, j)
        r += 1
        j -= 1
        yield r, (i, j)

【讨论】：

如果我理解正确的话，这是一条围绕 i,j 的方形蛇，对吗？只是警告杰森，以防他真的需要一个真正的半径。
赞成，因为它最接近我的要求，但它是正方形，而不是圆形螺旋。
@JasonSundram - 很好的编辑，明确地阐明了你在寻找什么。如果您可以手动生成您想要在圆圈周围最多几个循环的点并将其添加到您的问题中，那将非常有帮助，这将使尝试编写符合您期望的解决方案变得更加容易。

【解决方案5】：

虽然我不完全确定您要做什么，但我会这样开始：

def xxx():
    for row in M[i-R:i+R+1]:
        for val in row[j-R:j+r+1]:
            yield val

我不确定您希望螺旋的顺序是多少，这很重要吗？它是否必须按 R 顺序递增？还是从特定方位角开始顺时针？

R，曼哈顿的距离度量是多少？欧几里得？还有什么？

【讨论】：

【解决方案6】：

我要做的是使用阿基米德螺旋方程：

r(theta) = a + b*theta

然后将极坐标 (r,theta) 转换为 (x,y)，通过使用

x = r*cos(theta)
y = r*sin(theta)

cos 和 sin 在 math 库中。然后将得到的 x 和 y 舍入为整数。您可以将 x 和 y 偏移起始索引，以获得数组的最终索引。

但是，如果您只是对找到 f 返回 true 的第一个半径感兴趣，我认为执行以下伪代码会更有益：

for (i,j) in matrix:
    radius = sqrt( (i-i0)^2 + (j-j0)^2) // (i0,j0) is the "center" of your spiral
    radiuslist.append([radius, (i,j)])
sort(radiuslist) // sort by the first entry in each element, which is the radius
// This will give you a list of each element of the array, sorted by the
// "distance" from it to (i0,j0)
for (rad,indices) in enumerate(radiuslist):
    if f(matrix[indices]):
        // you found the first one, do whatever you want

【讨论】：

这确实会返回一个圆形螺旋，但它不会通过起点周围的所有网格点。我已经用一张图片更新了这个问题。

【解决方案7】：

好吧，我很尴尬，这是迄今为止我想出的最好的。但也许它会帮助你。由于它实际上不是循环迭代器，因此我不得不接受您的测试函数作为参数。

问题：

未针对跳过数组外的点进行优化
仍然使用方形迭代器，但它确实找到了最近的点
我没有使用过 numpy，所以它是为列表制作的。您需要更改的两点已注释
我将方形迭代器保留为长格式，以便于阅读。它可能更干燥

这里是代码。您问题的关键解决方案是顶级“spiral_search”函数，它在方形螺旋迭代器之上添加了一些额外的逻辑，以确保找到最近的点。

from math import sqrt

#constants
X = 0
Y = 1

def spiral_search(array, focus, test):
    """
    Search for the closest point to focus that satisfies test.
    test interface: test(point, focus, array)
    points structure: [x,y] (list, not tuple)
    returns tuple of best point [x,y] and the euclidean distance from focus
    """
    #stop if focus not in array
    if not _point_is_in_array(focus, array): raise IndexError("Focus must be within the array.")
    #starting closest radius and best point
    stop_radius = None
    best_point = None 
    for point in _square_spiral(array, focus):
        #cheap stop condition: when current point is outside the stop radius
        #(don't calculate outside axis where more expensive)
        if (stop_radius) and (point[Y] == 0) and (abs(point[X] - focus[X]) >= stop_radius):
            break #current best point is already as good or better so done
        #otherwise continue testing for closer solutions
        if test(point, focus, array):
            distance = _distance(focus, point)
            if (stop_radius == None) or (distance < stop_radius):
                stop_radius = distance
                best_point = point
    return best_point, stop_radius

def _square_spiral(array, focus):
    yield focus
    size = len(array) * len(array[0]) #doesn't work for numpy
    count = 1
    r_square = 0
    offset = [0,0]
    rotation = 'clockwise'
    while count < size:
        r_square += 1
        #left
        dimension = X
        direction = -1
        for point in _travel_dimension(array, focus, offset, dimension, direction, r_square):
            yield point
            count += 1
        #up
        dimension = Y
        direction = 1
        for point in _travel_dimension(array, focus, offset, dimension, direction, r_square):
            yield point
            count += 1
        #right
        dimension = X
        direction = 1
        for point in _travel_dimension(array, focus, offset, dimension, direction, r_square):
            yield point
            count += 1
        #down
        dimension = Y
        direction = -1
        for point in _travel_dimension(array, focus, offset, dimension, direction, r_square):
            yield point
            count += 1

def _travel_dimension(array, focus, offset, dimension, direction, r_square):
    for value in range(offset[dimension] + direction, direction*(1+r_square), direction):
        offset[dimension] = value
        point = _offset_to_point(offset, focus)
        if _point_is_in_array(point, array):
            yield point

def _distance(focus, point):
    x2 = (point[X] - focus[X])**2
    y2 = (point[Y] - focus[Y])**2
    return sqrt(x2 + y2)

def _offset_to_point(offset, focus):
    return [offset[X] + focus[X], offset[Y] + focus[Y]]

def _point_is_in_array(point, array):
    if (0 <= point[X] < len(array)) and (0 <= point[Y] < len(array[0])): #doesn't work for numpy
        return True
    else:
        return False

【讨论】：

我正在尝试查看并可视化人们为circle_around 返回的索引——有什么方法可以将此解决方案转换为只返回那些螺旋索引？
我不完全满意的原因是这并没有按照您的要求返回螺旋迭代。它使用平方迭代，但足够聪明，可以知道何时找到了实际最近的点（而不仅仅是第一个可能不是最近的点）。因此，如果您通过测试功能，它会为您找到第一点。我为其他几种方法编写了伪代码，但它们在计算沿实际圆或螺旋的 unique 点方面都很昂贵。传入你的测试函数会不会不够好？
我只是希望能够可视化您的结果，以便将它们与其他人的答案进行比较。我想我必须想出另一种方法来做到这一点。
如果你有自己的可视化方式，你只需要遍历_square_spiral()。