La forma más rápida de obtener todos los puntos entre dos coordenadas (X, Y) en python

Así que tengo un LineString shapely LineString :

 print np.round(shapely_intersecting_lines.coords).astype(np.int) >>> array([[ 1520, -1140], [ 1412, -973]]) 

Esto se puede interpretar como una matriz numpy como se ve arriba.

Quiero obtener todos los puntos intermedios, es decir, obtener los puntos de la línea intermedios como valores enteros. La salida debería ser algo como esto:

 array([[ 1520, -1140], [ 1519, -1139], [ 1519, -1138], ..., [ 1413, -975], [ 1412, -974], [ 1412, -973]], dtype=int32) 

Publiqué esto anteriormente en gis.stackexchange con la esperanza de que hubiera una solución shapely que fuera eficiente. La solución fue buena al principio, sin embargo, la solución ahora es demasiado lenta, ya que la ejecuté más de 50000 veces en mi código. En mi computadora, cada bucle toma aproximadamente 0.03 s, lo que resulta en más de un día de funcionamiento. Es demasiado lento para lo que necesito aquí y esperaba ver si alguien conoce una solución vectorizada para esto.

Bresenham puede ser inteligente, pero estoy bastante seguro de que la vectorización de la fuerza bruta es más rápida. He escrito dos variantes: la primera es más fácil de leer, la segunda es más rápida (80 us vs 50 us).

Actualización Se corrigió un error (gracias @Varlor) y se agregó una variante nd.

 import numpy as np from timeit import timeit def connect(ends): d0, d1 = np.abs(np.diff(ends, axis=0))[0] if d0 > d1: return np.c_[np.linspace(ends[0, 0], ends[1, 0], d0+1, dtype=np.int32), np.round(np.linspace(ends[0, 1], ends[1, 1], d0+1)) .astype(np.int32)] else: return np.c_[np.round(np.linspace(ends[0, 0], ends[1, 0], d1+1)) .astype(np.int32), np.linspace(ends[0, 1], ends[1, 1], d1+1, dtype=np.int32)] def connect2(ends): d0, d1 = np.diff(ends, axis=0)[0] if np.abs(d0) > np.abs(d1): return np.c_[np.arange(ends[0, 0], ends[1,0] + np.sign(d0), np.sign(d0), dtype=np.int32), np.arange(ends[0, 1] * np.abs(d0) + np.abs(d0)//2, ends[0, 1] * np.abs(d0) + np.abs(d0)//2 + (np.abs(d0)+1) * d1, d1, dtype=np.int32) // np.abs(d0)] else: return np.c_[np.arange(ends[0, 0] * np.abs(d1) + np.abs(d1)//2, ends[0, 0] * np.abs(d1) + np.abs(d1)//2 + (np.abs(d1)+1) * d0, d0, dtype=np.int32) // np.abs(d1), np.arange(ends[0, 1], ends[1,1] + np.sign(d1), np.sign(d1), dtype=np.int32)] def connect_nd(ends): d = np.diff(ends, axis=0)[0] j = np.argmax(np.abs(d)) D = d[j] aD = np.abs(D) return ends[0] + (np.outer(np.arange(aD + 1), d) + (aD>>1)) // aD ends = np.array([[ 1520, -1140], [ 1412, -73]]) ends_4d = np.array([[ 100, -302, 101, -49], [ -100, -45, 112, 100]]) print(connect(ends)) print(connect_nd(ends_4d)) assert np.all(connect(ends)==connect2(ends)) assert np.all(connect(ends)==connect_nd(ends)) assert np.all(connect(ends)==connect(ends[:, ::-1])[:, ::-1]) assert np.all(connect(ends)==connect(ends[::-1])[::-1]) print(timeit('f(ends)', globals={'f': connect, 'ends': ends}, number=10000)*100, 'us') print(timeit('f(ends)', globals={'f': connect2, 'ends': ends}, number=10000)*100, 'us') print(timeit('f(ends)', globals={'f': connect_nd, 'ends': ends}, number=10000)*100, 'us') 

Salida de muestra:

 [[ 1520 -1140] [ 1520 -1139] [ 1520 -1138] ..., [ 1412 -75] [ 1412 -74] [ 1412 -73]] [[ 100 -302 101 -49] [ 99 -301 101 -48] [ 98 -300 101 -48] ..., [ -98 -47 112 99] [ -99 -46 112 99] [-100 -45 112 100]] 78.8237597000034 us 48.02509490000375 us 62.78072760001123 us 

Usando generadores para que así ahorres memoria.

 import numpy as np import math d = np.array( [ [1520,-1140], [1412,-973] ],dtype=float); rise = d[0,1]-d[1,1]; run = d[0,0]-d[1,0]; slope = rise/run; print slope fromPt = d[0]; sign = 1 if (slope<0): sign = -1; points = ([d[0,0]+sign*i,math.floor(d[0,1]+(sign*i*slope))] for i in xrange(1+int(math.ceil(abs(d[0,0]-d[1,0]))))) for pt in points: print pt 

Aquí está el algoritmo de línea de Bresenhams como un generador. Si desea una lista, simplemente llame a list () en la salida:

 def line(x0, y0, x1, y1): deltax = x1-x0 dxsign = int(abs(deltax)/deltax) deltay = y1-y0 dysign = int(abs(deltay)/deltay) deltaerr = abs(deltay/deltax) error = 0 y = y0 for x in range(x0, x1, dxsign): yield x, y error = error + deltaerr while error >= 0.5: y += dysign error -= 1 yield x1, y1