python - lenguaje - matriz en espiral pseint

Iterar más de 2d matriz en una espiral circular en expansión

Question

Dada una n por n matriz M , en la fila i y columna j , me gustaría iterar sobre todos los valores vecinos en una espiral circular.

El punto de hacer esto es probar alguna función, f , que depende de M, para encontrar el radio lejos de (i, j) en el que f devuelve True . Entonces, f ve así:

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

y se llamaría así:

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

Donde circle_around es la función que devuelve (un iterador a) índices en una espiral circular. Entonces, para cada punto en M , este código calcularía y almacenaría el radio desde el punto en el que f devuelve True .

Si hay una forma más eficiente de calcular R , estaría abierto a eso también.

Actualizar:

Gracias a todos los que enviaron sus respuestas. He escrito una breve función para trazar el resultado de los iteradores circle_around , para mostrar lo que hacen. Si actualiza su respuesta o publica una nueva, puede usar este código para validar su solución.

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

Aquí están los resultados: plot(circle_around(0, 0), "FJ") :

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

He codificado la sugerencia de Magnesium de la siguiente manera:

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

Como puede ver, ninguno de los resultados que satisfacen la interfaz que busco produjo una espiral circular que cubre todos los índices alrededor de 0, 0. FJ es el más cercano, aunque WolframH golpea los puntos correctos, simplemente no en espiral orden.

Answer 1

Aquí hay una implementación basada en bucle para 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

Aunque no estoy del todo seguro de lo que estás tratando de hacer, comenzaría así:

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

No estoy seguro de cuánto orden quieres para tu espiral, ¿es importante? tiene que ser en orden R creciente? o tal vez en el sentido de las agujas del reloj a partir de azimut particular?

¿Cuál es la medida de distancia para R, manhattan? euclidiano? ¿algo más?

Answer 3

Bueno, estoy bastante avergonzado, este es el mejor que he encontrado hasta ahora. Pero tal vez te ayude. Como en realidad no es un iterador circular, tuve que aceptar su función de prueba como argumento.

Problemas:

• no optimizado para saltar puntos fuera de la matriz
• todavía usa un iterador cuadrado, pero encuentra el punto más cercano
• No he utilizado numpy, por lo que está hecho para la lista de listas. los dos puntos que necesita cambiar son comentados
• Dejé el iterador cuadrado en una forma larga para que sea más fácil de leer. podría ser más seco

Aquí está el código. La solución clave para su pregunta es la función de nivel superior "spiral_search" que agrega un poco de lógica adicional sobre el iterador espiral cuadrado para asegurarse de que se encuentre el punto más cercano.

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 4

Como se mencionó que el orden de los puntos no importa, simplemente los ordené por el ángulo ( arctan2 ) en el que aparecen en un radio dado. Cambia N para obtener más puntos.

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

Da para N=8 :

Más puntos N=16 (lo siento por el daltónico):

Esto se acerca claramente a un círculo y golpea cada punto de la cuadrícula en orden de radio creciente.

Answer 5

Lo que haría es usar la ecuación para una espiral de Arquímedes:

r(theta) = a + b*theta

y luego convierte las coordenadas polares (r, theta) en (x, y), usando

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

cos y el sin están en la biblioteca de math . Luego redondee los números x e y resultantes a enteros. Puede compensar xey después con el índice inicial, para obtener los índices finales de la matriz.

Sin embargo, si solo está interesado en encontrar el primer radio donde f devuelve verdadero, creo que sería más beneficioso hacer el siguiente pseudocódigo:

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

Si sigue los índices helicoidales xey, observará que ambos pueden definirse de forma recursiva. Por lo tanto, es bastante fácil crear una función que genere recursivamente los índices correctos:

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

Como puede ver en la imagen de arriba, esto proporciona exactamente lo que se solicita.

Answer 7

Una manera de ceder puntos al aumentar la distancia es dividirla en partes fáciles y luego combinar los resultados de las partes. Es bastante obvio que itertools.merge debería hacer la fusión. Las partes fáciles son columnas , porque para x fijo, los puntos (x, y) se pueden ordenar mirando solo el valor de y.

A continuación se muestra una implementación (simplista) de ese algoritmo. Tenga en cuenta que se utiliza la distancia euclidiana al cuadrado y que se incluye el punto central. Lo que es más importante, solo se consideran los puntos (x, y) con x en el range(x_end) , pero creo que está bien para su caso de uso (donde x_end sería n en su notación anterior).

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

Editar: Código de prueba:

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

Resultado de la prueba (los puntos están numerados en el orden en que aparecen en 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                

Edición 2: si realmente necesita valores negativos de i , reemplace el range(end_x) con el range(-end_x, end_x) en la función cirlce_around .