Итерация по 2d массиву в расширяющейся круговой спирали

Question

Учитывая n на 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"):

Я закодировал предложение Магния следующим образом:

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 достигает правильных точек , только не по спирали.

Answer 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 (извините за дальтоник):

Это ясно приближается к кругу и поражает каждую точку сетки в порядке увеличения радиуса.

Answer 2

Один из способов получения точек с увеличением расстояния - разбить его на легких частей , а затем объединить результатов частей вместе. Совершенно очевидно, что itertools.merge должно выполнить слияние. easy parts - это столбцы , поскольку для фиксированного x точки (x, y) можно упорядочить, взглянув только на значение y.

Ниже приведена (упрощенная) реализация этого алгоритма. Обратите внимание, что используется евклидово расстояние в квадрате и указывается центральная точка. Что наиболее важно, учитываются только точки (x, y) с x в range(x_end), но я думаю, что это нормально для вашего варианта использования (где 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, замените range(end_x) на range(-end_x, end_x) в функции cirlce_around.

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

Как видно из изображения выше, это дает именно то, чтоспросил.

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

Answer 5

Я бы использовал уравнение для архимедовой спирали:

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

Answer 6

Хотя я не совсем уверен, что вы пытаетесь сделать, я бы начал так:

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

Я не уверен, какой порядок вы хотите для своей спирали, это вообще важно? это должно быть в порядке увеличения R? или, возможно, по часовой стрелке, начиная с определенного азимута?

Что такое мера расстояния для R, Манхэттен? евклидовым? что-то еще?

Answer 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