Penetração de sinal numpy

O problema consiste em modelar a propagação de um sinal (por exemplo, luz ou som etc.) através de uma série de obstáculos, como na figura abaixo. O sinal não pode passar pela superfície inferior (terreno), mas pode atravessar obstáculos. Eu quero contar o número de obstáculos atravessados.

Terreno e obstáculos estão em matrizes numpy 2D (x, y, z). Isto é o que eu faço:

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

O resultado seria algo como [0, 0, 0, 1, 1, 1, 2, 3, 4, 4, 4 ...]por linha.

Este método funciona bem (desde que os vales no terreno sejam preenchidos usando numpy.maximum.accumulate). Agora, seria possível acelerar a coisa usando uma solução vetorizada?

Como os comentários dizem acima, você provavelmente pode vetorizar a operação para remover os loops for e torná-la mais eficiente.

No entanto, se você considerar o problema de uma maneira um pouco diferente - a de limiar - poderá tirar proveito das ferramentas do scipy ndimage para contar os obstáculos:

Primeiro, limite os dados do terreno pela altura do sinal para obter uma matriz booleana de onde o sinal poderia estar, independentemente da origem.

signal_reach = terrain < signal_height

Em seguida, você pode usar o ndimage.labelmétodo para agrupar regiões distintas:

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

Uma vez feito isso, obtenha o ID da região que corresponde à (s) célula (s) de origem do seu sinal. No seu caso, seria a primeira coluna.

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]

Agora, o limite é onde o sinal intercepta o terreno + obstáculos e, desta vez, filtra regiões fora ou áreas de interesse usando numpy.isin.

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

E uma rodada final ndimage.labelmostra os obstáculos interceptados, porque já filtramos as áreas bloqueadas pelo terreno:

obstacles_hit, obstacle_count = ndimage.label(signal_intercepted)

Há um pouco mais de código aqui, mas há duas grandes vantagens:

• Não para loops significa que o código deve ser razoavelmente rápido,
• E, para fins de cálculo, a origem do sinal pode estar em várias células, em qualquer lugar da varredura do terreno.