Señal de penetración en numpy

El problema consiste en modelar la propagación de una señal (por ejemplo, luz o sonido, etc.) a través de una serie de obstáculos, como en la figura a continuación. La señal no puede pasar a través de la superficie inferior (terreno), pero puede atravesar obstáculos. Quiero contar la cantidad de obstáculos atravesados.

El terreno y los obstáculos están en matrices numpy 2D (x, y, z). Esto es lo que hago:

output = numpy.zeros(terrain.shape)

obstacles = terrain + obstacle_heights


for i in xrange (obstacles.shape[0]):
    for j in xrange (obstacles.shape[1]):

        mask = obstacles[i,j] > terrain[i,j:]
        output[i,j:][mask] +=1

El resultado sería algo así [0, 0, 0, 1, 1, 1, 2, 3, 4, 4, 4 ...]por fila.

Este método funciona bien (siempre que los valles del terreno se llenen usando numpy.maximum.accumulate). Ahora, ¿sería posible acelerar la cosa usando una solución vectorizada?

Como dicen los comentarios anteriores, probablemente podría vectorizar la operación para eliminar los bucles for y hacerlo más eficiente.

Sin embargo, si considera el problema de una manera ligeramente diferente, la de los umbrales, puede aprovechar las herramientas de nipimage scipy para contar los obstáculos:

Primero, limite los datos del terreno por la altura de su señal para obtener una matriz booleana de dónde podría estar la señal, independientemente del origen.

signal_reach = terrain < signal_height

Luego puede usar el ndimage.labelmétodo para agrupar regiones discretas:

from scipy import ndimage
signal_regions, region_count = ndimage.label(signal_reach)

Una vez hecho esto, obtenga las ID de región que coinciden con las celdas de origen de su señal. En su caso, sería la primera columna.

import numpy as np
origin_labels = np.unique(signal_regions[:, 0]) # or whatever indexes meet the start of the signal
# ndimage lables are greater than 1, 0 is an unlabeled region
origin_labels = origin_lables[origin_lables > 0]

Ahora el umbral es donde la señal intercepta el terreno + obstáculos, y esta vez filtra las regiones externas o el área de interés utilizando numpy.isin.

area_of_interest = np.isin(signal_regions, origin_labels)
signal_intercepted = (obstacles >= signal_height) & area_of_interest

Y una ronda final de ndimage.labelle da una cuenta de los obstáculos interceptados, porque ya hemos filtrado las áreas bloqueadas por el terreno:

obstacles_hit, obstacle_count = ndimage.label(signal_intercepted)

Hay un poco más de código aquí, pero hay dos grandes ventajas:

• No para bucles significa que el código debe ser razonablemente rápido,
• Y para fines de cálculo, el origen de la señal puede estar en muchas celdas, en cualquier parte del ráster del terreno.