python 中的有界圆插值

Question

我的问题与这里的问题类似。简而言之，我有一个时间序列角度数据，其范围在 [0, 360] 之间。我需要计算测量之间的插值。目前，我正在使用 scipy.interpolate.interp1d。为了清楚地说明我的问题，这里有一个例子，

import numpy as np
from scipy import interpolate
data = np.array([[0, 2, 4], [1, 359, 1]]) # first row time index, second row angle measurements
f = interpolate.interp1d(data[0, :], data[1, :], kind='linear', bounds_error=False, fill_value=None)
f([1, 3])

这将导致 [ 180., 180.]。然而，在时间 2 和时间 4 之间，角度从 359 变为 1，这只是 2 度的变化，并且 3 处的插值应为 0。随着时间的推移，角度在 CCW 方向上发生变化。

最后，我的问题是，

我可以使用任何标准模块来实现此目的吗？

只是因为我想尽可能避免自定义方法！

Answer 1

每次检测到存在跳跃时，只需添加 360° 补数，然后使用模运算恢复到第一个 360 度。例如：

In [1]: import numpy as np

In [2]: from scipy import interpolate

In [3]: data = np.array([[0, 2, 4, 6, 8], [1, 179, 211, 359, 1]])

In [4]: complement360 = np.rad2deg(np.unwrap(np.deg2rad(data[1])))

In [5]: complement360
Out[5]: array([   1.,  179.,  211.,  359.,  361.])

In [6]: f = interpolate.interp1d(data[0], complement360, kind='linear', bounds_error=False, fill_value=None)

In [7]: f(np.arange(9))
Out[7]: array([   1.,   90.,  179.,  195.,  211.,  285.,  359.,  360.,  361.])

In [8]: f(np.arange(9))%360
Out[8]: array([   1.,   90.,  179.,  195.,  211.,  285.,  359.,    0.,    1.])

备注，我确实在这里添加了一些额外的值，否则

np.unwrap

没有现实的方法可以知道角度在哪个方向增加，这也可能是你知道角度以这种方式增加的方式（差异连续值之间的角度小于 180°，除非存在实际的不连续性）。

但是，如果您确实拥有在两个连续项目之间角度跳跃大于 180° 的数据，但您知道角度变化的方向（例如 CCW）并且它是单调变化的，那么您可以像这样检测它：

In [31]: data = np.array([1, 359, 1, 60, 359, 177, 2])  # mock-data

In [32]: jumps = np.diff(data)<0  # assumptions: angle increases stricly monotonously CCW

In [33]: np.hstack((data[0], data[1:] + np.cumsum(np.sign(d)<0)*360))
Out[33]: array([   1,  359,  361,  420,  719,  897, 1082])

Answer 2

从版本 1.10.0 开始，numpy.interp 采用句点关键字：http://docs.scipy.org/doc/numpy/reference/ generated/numpy.interp.html

Answer 3

我认为下面的代码完全符合您的要求。首先，它展开角度，使得连续值之间的跳跃永远不会大于 180 度（请参阅 MATLAB 关于展开的精彩文档），然后执行插值，最后使用

np.mod

使得返回的角度位于区间 [0, 360)。

import numpy as np

def circular_interpolation(x: np.ndarray, xp: np.ndarray, fp: np.ndarray, period: float = 360) -> np.ndarray:
    """
    One dimensional linear interpolation for monotonically increasing sample points where points first are unwrapped,
    secondly interpolated and finally bounded within the specified period.
    
    :param x: The x-coordinates at which to evaluate the interpolated values.
    :param xp: The x-coordinates of the data points, must be increasing.
    :param fp: The y-coordinates of the data points, same length as `xp`.
    :param period: Size of the range over which the input wraps.
    :return: The interpolated values, same shape as `x`.
    """
    
    y = np.mod(np.interp(x, xp, np.unwrap(fp, period=period)), period)
    return y

---

请注意，

period

的

np.unwrap

选项是 Numpy 1.21.0 版本中的新选项。对于旧版本的 Numpy，您可以在执行插值之前转换为弧度。

Answer 4

scipy 提供了 Slerp，可用于插值 3D 旋转。角度可以视为一维旋转，是三维旋转的特例。所以我们可以利用scipy的旋转插值来做角度插值。下面是一个例子：

import numpy as np
from scipy.spatial.transform import Rotation as R
from scipy.spatial.transform import Slerp`    

def slerp_angles(angle1, angle2, num_points):
    rot1 = R.from_euler('xyz', [angle1, 0, 0])
    rot2 = R.from_euler('xyz', [angle2, 0, 0])
    key_times = [0, 1]
    rotations = R.concatenate([rot1, rot2])
    slerp = Slerp(key_times, rotations)
    interp_times = np.linspace(0, 1, num_points)
    interp_rotations = slerp(interp_times)
    angles = interp_rotations.as_euler('xyz')[:, 0]
    return angles

num_points = 10

angle1, angle2 = (- np.pi + 0.01, np.pi - 0.01)
interpolated_angles = slerp_angles(angle1, angle2, num_points)
print(f"Interpolating from {angle1} to {angle2}:")
print(interpolated_angles)

输出：

Interpolating from -3.1315926535897933 to 3.1315926535897933:
[-3.13159265 -3.13381488 -3.1360371  -3.13825932 -3.14048154  3.14048154
  3.13825932  3.1360371   3.13381488  3.13159265]