Iterate over 2d array in an expanding circular spiral

Question

Given an n by n matrix M, at row i and column j, I'd like to iterate over all the neighboring values in a circular spiral.

The point of doing this is to test some function, f, which depends on M, to find the radius away from (i, j) in which f returns True. So, f looks like this:

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

and would be called like this:

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

Where circle_around is the function that returns (an iterator to) indices in a circular spiral. So for every point in M, this code would compute and store the radius from that point in which f returns True.

If there's a more efficient way of computing R, I'd be open to that, too.

---

Update:

Thanks to everyone who submitted answers. I've written a short function to plot the output from your circle_around iterators, to show what they do. If you update your answer or post a new one, you can use this code to validate your solution.

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")

Here are the results: plot(circle_around(0, 0), "F.J"):

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

I've coded up Magnesium's suggestion as follows:

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"):

As you can see, none of the results that satisfy the interface I'm looking for have produced a circular spiral that covers all of the indices around 0, 0. F.J's is the closest, although WolframH's hits the right points, just not in spiral order.

Answer 1

Here is a loop based implementation for 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)

Answer 2

If you follow the x and y helical indices you notice that both of them can be defined in a recursive manner. Therefore, it is quite easy to come up with a function that recursively generates the correct indices:

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")

As you can see from the image above, this gives exactly what's asked for.

Answer 3

One way for yielding points with increasing distance is to break it down into easy parts, and then merge the results of the parts together. It's rather obvious that itertools.merge should do the merging. The easy parts are columns, because for fixed x the points (x, y) can be ordered by looking at the value of y only.

Below is a (simplistic) implementation of that algorithm. Note that the squared Euclidian distance is used, and that the center point is included. Most importantly, only points (x, y) with x in range(x_end) are considered, but I think that's OK for your use case (where x_end would be n in your notation above).

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

Edit: Test code:

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)))

Result of test (points are numbered in the order they are yielded by 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                

Edit 2: If you really do need negative values of i, replace range(end_x) with range(-end_x, end_x) in the cirlce_around function.

Answer 4

What I would do is use the equation for an Archimedean spiral:

r(theta) = a + b*theta

and then convert the polar coordinates (r,theta) into (x,y), by using

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

cos and sin are in the math library. Then round the resulting x and y to integers. You can offset x and y afterward by the starting index, to get the final indices of the array.

However, if you are just interested in finding the first radius where f returns true, I think it would be more beneficial to do the following pseudocode:

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

Answer 5

Well, I'm pretty embarrassed this is the best I have come up with so far. But maybe it will help you. Since it's not actually a circular iterator, I had to accept your test function as an argument.

Problems:

• not optimized to skip points outside the array
• still uses a square iterator, but it does find the closest point
• i haven't used numpy, so it's made for list of lists. the two points you need to change are commented
• i left the square iterator in a long form so it's easier to read. it could be more DRY

Here is the code. The key solution to your question is the top level "spiral_search" function which adds some extra logic on top of the square spiral iterator to make sure that the closest point is found.

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

Answer 6

Although I'm not entirely sure what you are trying to do, I'd start like this:

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

I'm not sure how much order you want for your spiral, is it important at all? does it have to be in increasing R order? or perhaps clockwise starting at particular azimuth?

What is the distance measure for R, manhattan? euclidean? something else?

Answer 7

Since it was mentioned that the order of the points do not matter, I've simply ordered them by the angle (arctan2) in which they appear at a given radius. Change N to get more points.

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()

Gives for N=8:

More points N=16 (sorry for the colorblind):

This clearly approaches a circle and hits every grid point in order of increasing radius.