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.

Solution 1

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.

Solution 2

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.

Solution 3

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.

Solution 4

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)

View more solutions

Comments

[KobeJohn]

Can you confirm that your arrays are very large or you have to do this many times or the truth function is expensive or...? I could come up with something, but it seems like premature optimization unless you really need to avoid testing outside the radius. Of course the simple solution would be to find the radius for every false point in the array and then just find the smallest radius. Neat problem though if you really need it.

[pythonBOI]

@yakiimo, the arrays have 1-2 million entries.

[KobeJohn]

Does F.J's answer work for you or do you need a real circle?

[pythonBOI]

I would prefer a real circle.

[KobeJohn]

FYI I still have this on my to-do list but it's quite long at the moment.

[pythonBOI]

@yakiimo -- awesome. I've been meaning to write some code to visualize the generated indices to show their spiral-ness (or lack thereof), but I haven't gotten around to it yet.

[KobeJohn]

similar (basically same) question but still not calculating based on a circular radius

[KobeJohn]

I believe the code I posted will give you the answer you want although in a slightly different way than you asked for. If using a square iterator, it's really important to handle the case of the first point not necessarily being the closest due to the square shape of the iterator.

[Hooked]

Is the ordering of the points important in the spiral, or do you simply want all points (i,j) such that f(i,j) < r?

[pythonBOI]

@Hooked, I want them ordered by ascending radius. Within any set of points at the same distance, I suppose I don't care what order they are returned in, although it would be nice if it was a nice order.

[Reinstate Monica]

@JasonSundram: Why does your matrix have negative indices?

[John]

Related math.stackexchange question: math.stackexchange.com/questions/1740130/…

[KobeJohn]

If I understand correctly, this is a square snake around i,j correct? Just a warning to Jason in case he actually needs a real radius.

[pythonBOI]

I'm trying to go through and visualize the indices people are returning for circle_around -- is there any way you can convert this solution to just return those spiral indices?

[pythonBOI]

upvoted since it's the closest to what I've asked for, but it is a square, not a circular spiral.

[pythonBOI]

this does return a circular spiral, but it doesn't pass through all the grid points around the starting point. I've updated the question with a picture of what this produces.

[pythonBOI]

this doesn't look anything like a spiral -- see my updates to the question. did I miss something? What got produced was this: i.stack.imgur.com/h0mNa.png

[Andrew Clark]

@JasonSundram - Nice edit, definitely clarifies what you're looking for. It would be really helpful if you could manually generate the points you would want up to a few loops around the circle and add that to your question as well, that would make it a lot easier to try to code a solution that matches your expectation.

[Reinstate Monica]

@JasonSundram: 1) I assumed nonnegative indices, because you wrote about an n x n-Matrix. 2) I also assumed that the order of points of the same distance is not important. That's the impression I got from your question and comments. Is that wrong?

[pythonBOI]

gotcha. Thanks for the clarification and the code to show how your answer works. I've also regenerated the image for your code in the post to more accurately depict how it works.

[KobeJohn]

The reason I wasn't totally satisfied is that this doesn't return a spiral iteration as you asked for. It uses a square iteration but is smart enough to know when the actual closest point has been found (rather than just the first which may not be the closest). So it will work for you to find the first point if you pass it your test function. I did pseudocode for several other ways but they were all expensive in terms of calculating unique points that are along an actual circle or spiral. Will passing in your test function not work well enough?

[pythonBOI]

I just wanted to be able to visualize your results to compare them to everyone else's answers. I guess I'll have to come up with another way to do that.

[KobeJohn]

if you have your own way to visualize, all you need is to iterate through _square_spiral().