Оптимизация евклидовой функции Numpy для расстояния и направления

Question

Я пытаюсь вычислить евклидово расстояние и направление из исходной координаты в массиве.

Графический пример

Вот что я былв состоянии придумать, однако это довольно медленно для больших массивов.Евклидово расстояние и направление, основанные на исходных координатах, сильно зависят от индекса каждой ячейки.Вот почему я зацикливаю каждую строку и столбец.Я посмотрел на scipy cdist, pdist и np linalg.

import numpy as np
from math import atan, degrees, sqrt
from timeit import default_timer

def euclidean_from_source(input_array, y_index, x_index):
    # copy arrays
    distance = np.empty_like(input_array, dtype=float)
    direction = np.empty_like(input_array, dtype=int)

    # loop each row
    for i, row in enumerate(X):
        # loop each cell
        for c, cell in enumerate(row):
            # get b
            b = x_index - i
            # get a
            a = y_index - c

            hypotenuse = sqrt(a * a + b * b) * 10
            distance[i][c] = hypotenuse
            direction[i][c] = get_angle(a, b)

    return [distance, direction]

def calibrate_angle(a, b, angle):
    if b > 0 and a > 0:
        angle+=90
    elif b < 0 and a < 0:
        angle+=270
    elif b > 0 > a:
        angle+=270
    elif a > 0 > b:
        angle+=90
    return angle

def get_angle(a, b):
    # get angle
    if b == 0 and a == 0:
        angle = 0
    elif b == 0 and a >= 0:
        angle = 90
    elif b == 0 and a < 0:
        angle = 270
    elif a == 0 and b >= 0:
        angle = 180
    elif a == 0 and b < 0:
        angle = 360
    else:
        theta = atan(b / a)
        angle = degrees(theta)

    return calibrate_angle(a, b, angle)

if __name__ == "__main__":
    dimension_1 = 5
    dimension_2 = 5

    X = np.random.rand(dimension_1, dimension_2)
    y_index = int(dimension_1/2)
    x_index = int(dimension_2/2)

    start = default_timer()
    distance, direction = euclidean_from_source(X, y_index, x_index)
    print('Total Seconds {{'.format(default_timer() - start))

    print(distance)
    print(direction)

ОБНОВЛЕНИЕ Я смог использовать функцию трансляции, чтобы сделать именно то, что мне нужно, и на долю отскорость.однако я все еще выясняю, как откалибровать угол к 0, 360 по всей матрице (модуль не будет работать в этом сценарии).

import numpy as np
from math import atan, degrees, sqrt
from timeit import default_timer


def euclidean_from_source_update(input_array, y_index, x_index):
    size = input_array.shape
    center = (y_index, x_index)

    x = np.arange(size[0])
    y = np.arange(size[1])

    # use broadcasting to get euclidean distance from source point
    distance = np.multiply(np.sqrt((x - center[0]) ** 2 + (y[:, None] - center[1]) ** 2), 10)

    # use broadcasting to get euclidean direction from source point
    direction = np.rad2deg(np.arctan2((x - center[0]) , (y[:, None] - center[1])))

    return [distance, direction]

def euclidean_from_source(input_array, y_index, x_index):
    # copy arrays
    distance = np.empty_like(input_array, dtype=float)
    direction = np.empty_like(input_array, dtype=int)

    # loop each row
    for i, row in enumerate(X):
        # loop each cell
        for c, cell in enumerate(row):
            # get b
            b = x_index - i
            # get a
            a = y_index - c

            hypotenuse = sqrt(a * a + b * b) * 10
            distance[i][c] = hypotenuse
            direction[i][c] = get_angle(a, b)
    return [distance, direction]

def calibrate_angle(a, b, angle):
    if b > 0 and a > 0:
        angle+=90
    elif b < 0 and a < 0:
        angle+=270
    elif b > 0 > a:
        angle+=270
    elif a > 0 > b:
        angle+=90
    return angle

def get_angle(a, b):
    # get angle
    if b == 0 and a == 0:
        angle = 0
    elif b == 0 and a >= 0:
        angle = 90
    elif b == 0 and a < 0:
        angle = 270
    elif a == 0 and b >= 0:
        angle = 180
    elif a == 0 and b < 0:
        angle = 360
    else:
        theta = atan(b / a)
        angle = degrees(theta)

    return calibrate_angle(a, b, angle)

if __name__ == "__main__":
    dimension_1 = 5
    dimension_2 = 5

    X = np.random.rand(dimension_1, dimension_2)
    y_index = int(dimension_1/2)
    x_index = int(dimension_2/2)

    start = default_timer()
    distance, direction = euclidean_from_source(X, y_index, x_index)
    print('Total Seconds {}'.format(default_timer() - start))

    start = default_timer()
    distance2, direction2 = euclidean_from_source_update(X, y_index, x_index)
    print('Total Seconds {}'.format(default_timer() - start))

    print(distance)
    print(distance2)

    print(direction)
    print(direction2)

Обновление 2 Спасибо всем за ответыпосле методов тестирования эти два метода были самыми быстрыми и дали нужные мне результаты.Я по-прежнему открыт для любых оптимизаций, которые вы, ребята, можете придумать.

def get_euclidean_direction(input_array, y_index, x_index):
    rdist = np.arange(input_array.shape[0]).reshape(-1, 1) - x_index
    cdist = np.arange(input_array.shape[1]).reshape(1, -1) - y_index
    direction = np.mod(np.degrees(np.arctan2(rdist, cdist)), 270)

    direction[y_index:, :x_index]+= -90
    direction[y_index:, x_index:]+= 270
    direction[y_index][x_index] = 0

    return direction

def get_euclidean_distance(input_array, y_index, x_index):
    size = input_array.shape
    center = (y_index, x_index)

    x = np.arange(size[0])
    y = np.arange(size[1])
    return np.multiply(np.sqrt((x - center[0]) ** 2 + (y[:, None] - center[1]) ** 2), 10)

Answer 1

Может быть проще просто передать коррдинат, который вы хотите измерить как массив или кортеж.Кроме того, хотя это может занять немного больше памяти, я думаю, что использование np.indices может быть немного быстрее для вычислений (поскольку оно позволяет np.einsum творить чудеса).

def euclidean_from_source(input_array, coord):    
    grid = np.indices(input_array.shape)
    grid -= np.asarray(coord)[:, None, None]
    distance = np.einsum('ijk, ijk -> jk', grid, grid) ** .5
    direction = np.degrees(np.arctan2(grid[0], grid[1]))
    return distance, direction

Этот методтакже немного более расширяемый для nd (хотя очевидно, что вычисления углов были бы немного сложнее

Answer 2

Эта операция чрезвычайно проста для векторизации.С одной стороны, a и b вообще не нужно вычислять в 2D, поскольку они зависят только от одного направления в массиве.Расстояние можно вычислить с помощью np.hypot.Вещание преобразует фигуру в правильную 2D-форму.

Ваша функция угла почти точно эквивалентна применению np.degrees к np.arctan2.

Непонятно, почему вы помечаете свои строки с помощью x, а столбцы с помощью y вместо стандартного способа сделать это, но до тех пор, пока вы последовательны, все должно быть в порядке.

Так вотвекторизованная версия:

def euclidean_from_source(input_array, c, r):
    rdist = np.arange(input_array.shape[0]).reshape(-1, 1) - r
    # Broadcasting doesn't require this second reshape
    cdist = np.arange(input_array.shape[1]).reshape(1, -1) - c
    distance = np.hypot(rdist, cdist) * 10
    direction = np.degrees(np.arctan2(rdist, cdist))
    return distance, direction

Я оставлю читателю в качестве упражнения определение того, требуется ли какая-либо дополнительная обработка для точной настройки угла, и, если да, то ее векторизация.