代码之家  ›  专栏  ›  技术社区  ›  dfresh22

优化欧氏距离和方向函数

  •  1
  • dfresh22  · 技术社区  · 8 年前

    我试图从numpy数组中的源坐标计算欧几里德距离和方向。

    图形示例 Example of the desired output

    这是我能够想到的,但是对于大型阵列来说它相对较慢。基于源坐标的欧氏距离和方向在很大程度上依赖于每个单元的索引。这就是为什么我要循环每一行和每一列。我研究了scipycdist、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)
    
    2 回复  |  直到 8 年前
        1
  •  2
  •   Daniel F    8 年前

    这个操作非常容易矢量化。首先, a b 根本不需要在2D中计算,因为它们只依赖于数组中的一个方向。距离可以用 np.hypot . 广播将把形状转换成正确的二维形式。

    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
    

    我将把它作为一个练习留给读者,以确定是否需要任何额外的处理来微调角度,如果需要,以矢量化的方式实现它。

        2
  •  1
  •   Daniel F    8 年前

    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
    

    这种方法对n-d的可扩展性也更强(尽管很明显角度计算会更复杂一些)