在 ndarray 的特定轴上向量化凸包和插值循环

Question

我正在努力寻找一种有效的方法来实现这种插值凸包数据处理。

我有一个 2D ndarray，称之为

arr

，形状为 (2000000,19)，包含浮点数。 我有一个 1D ndarray，称之为

w

，形状为 (19,)，也包含浮点数。

我所取得的成就（并且工作完美，除了速度极其缓慢）如下：

import numpy as np
from scipy.interpolate import interp1d

# Sample data
arr = np.array([[49.38639913, 49.76769437, 49.66370476, 49.49905455, 49.15242251,
        48.0518658 , 45.998071  , 45.31347273, 45.29614113, 45.25281212,
        45.0448329 , 44.61154286, 43.72763117, 42.38443203, 41.17121991,
        40.48662165, 40.35663463, 39.88001558, 39.55938095],
       [47.97387359, 47.86121818, 47.69656797, 47.18528571, 46.70000087,
        45.39146494, 43.50232035, 43.18168571, 43.82295498, 43.62364156,
        43.31167273, 42.88704848, 42.37576623, 41.0585645 , 40.37396623,
        39.09142771, 38.79679048, 38.51948485, 38.52815065]])
w = np.array([2.1017, 2.1197, 2.1374, 2.1548, 2.172 , 2.1893, 2.2068, 2.2254,
       2.2417, 2.2592, 2.2756, 2.2928, 2.3097, 2.326 , 2.3421, 2.3588,
       2.3745, 2.3903, 2.4064])

def upper_andrews_hull(points: np.ndarray):
    """
    Computes the upper half of the convex hull of a set of 2D points.
    :param points: an iterable sequence of (x, y) pairs representing the points.
    """
    # 2D cross product of OA and OB vectors, i.e. z-component of their 3D cross product.
    # Returns a positive value, if OAB makes a counter-clockwise turn,
    # negative for clockwise turn, and zero if the points are collinear.
    def cross(o, a, b):
        return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

    # Reverse the points so that we can pop from the end
    points = np.flip(points, axis=0)
    # Build upper hull
    upper = []
    for p in points:
        while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
            upper.pop()
        upper.append(p)
    # Reverse the upper hull
    upper = np.flip(np.array(upper), axis=0)

    return upper

result = np.array(arr.shape)
for i in range(arr.shape[0]):
    # Create points, using w as x values, and arr as y values
    points = np.stack((w, arr[i,:]), axis=1)
    # Calculate the convex hull around the points
    hull = upper_andrews_hull(points)
    # Interpolate the hull
    interp_function = interp1d(*hull.T)
    # Store interpolation's result to have the same x references as original points
    result[i,:] = interp_function(w)

我很确定有一种方法可以放弃循环而只使用向量计算，但我找不到它（另外，还有一个问题是

hull

并不总是具有相同数量的点，所以所有的船体无法存储在 ndarray 中。

Answer 1

这是找到上凸包的更快方法。

1. 使用

   spicy.spatial.ConvexHull

   找到凸包
2. 按逆时针顺序，从 x 最大的点到 x 最小的点取凸包的所有顶点

from scipy.spatial import ConvexHull

def upper_hull(points):
    ch = ConvexHull(points)
    max_idx = points[ch.vertices].argmax(axis=0)[0]
    min_idx = points[ch.vertices].argmin(axis=0)[0]
    ch_idxs = ch.vertices.tolist()
    upper_ch_idxs = (
        ch_idxs[max_idx:] + ch_idxs[:min_idx + 1] if min_idx < max_idx
        else ch_idxs[max_idx:min_idx + 1]
    )
    return points[upper_ch_idxs]