Optimizing Numpy Euclidean Distance and Direction Function

Question

I am trying to calculate the euclidean distance and direction from a source coordinate within a numpy array.

Graphic Example

Here is what I was able to come up with, however it is relatively slow for large arrays. Euclidean Distance and Direction based on source coordinates rely heavily on the index of each cell. that is why I am looping each row and column. I have looked into scipy cdist, pdist, and 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)

UPDATE I was able to use the broadcasting function to do exactly what I needed, and at a fraction of the speed. however I am still figuring out how to calibrate the angle to 0, 360 throughout the matrix (modulus will not work in this scenario).

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)

Update 2 Thanks everyone for the responses, after testing methods, these two methods were the fastest and produced the results I needed. I am still open to any optimizations you guys can think of.

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)

Are `x_index` and `y_index` always integers? Also, are there bounds on `input_array.shape`? You may be able to pre-comupte and use the same results for multiple runs. — Daniel F, Sep 21 '18 at 06:56
yup, they are always integers, and the maximum bounds I would set at 10000, 10000, what are you thinking? pre-compute and pickle the object? — dfresh22, Sep 21 '18 at 07:18
Possible duplicate of [How can the Euclidean distance be calculated with NumPy?](https://stackoverflow.com/questions/1401712/how-can-the-euclidean-distance-be-calculated-with-numpy) — user2699, Sep 21 '18 at 12:44
You'll want to select an answer to take your question off the unanswered queue at some point (assuming of course that you were able to solve your problem using the information provided). — Mad Physicist, Sep 23 '18 at 14:58

score 2 · Accepted Answer · edited Sep 21 '18 at 05:49

This operation is extremely easy to vectorize. For one thing, a and b don't need to be computed in 2D at all, since they only depend on one direction in the array. The distance can be computed with np.hypot. Broadcasting will convert the shape into the correct 2D form.

Your angle function is almost exactly equivalent to applying np.degrees to np.arctan2.

It's unclear why you label your rows with xand columns with y instead of the standard way of doing it, but as long as you're consistent, it should be fine.

So here's the vectorized version:

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

I will leave it as an exercise for the reader to determine if any additional processing is necessary to fine-tune the angle, and if so, to implement it in a vectorized manner.

looks like I updated the question at the same time you posted haha. I was using broadcasting, but most likely your method is faster. I will give it a shot — dfresh22, Sep 21 '18 at 05:29

score 1 · Answer 2 · answered Sep 21 '18 at 06:26

Might be easier to just pass the corrdinate you want to measure as an array or tuple. Also, while it might take a bit more memory, I think using np.indices might be a bit faster for the calculation (as it allows np.einsum to do its magic).

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

This method is also a bit more extensible to n-d (although obviously the angle calculations would be a bit trickier

I will give it a shot and compare the methods, so my broadcasting method was the fastest, but I will test yours. Thank you. I am using size (10000. 10000) to test — dfresh22, Sep 21 '18 at 06:30

Optimizing Numpy Euclidean Distance and Direction Function

2 Answers2