Python中的圆插值

Posted 2023-02-18

【中文标题】Python中的圆插值

【英文标题】：Circular interpolation in Python

【发布时间】：2013-03-31 02:45:36

【问题描述】：

我有两个系统，每个系统都有一个方向传感器（0-360 度），但是根据每个系统的方向和每个传感器的线性度，传感器可以提供截然不同的值。我有一个机械参考，可用于生成每个系统实际指向的表格。这会产生一个包含三列的表：

Physical  SystemA  SystemB
--------  -------  -------
 000.0     005.7    182.3
 005.0     009.8    178.4
 ...       ...      ...

仅从显示的数据中，我们可以看到 SystemA 离物理参考不远，但 SystemB 大约偏离 180 度，并且方向相反（想象一下它是倒置安装的）。

我需要能够在所有三个值之间来回映射：如果 SystemA 报告某物位于 105.7，我需要告诉用户那是什么物理方向，然后告诉 SystemB 指向同一个位置。如果 SystemB 进行初始报告，则相同。并且用户可以请求两个系统都指向所需的物理方向，因此需要告知 SystemA 和 SystemB 指向的位置。

线性插值并不难，但当数据方向相反时我遇到了麻烦，并且是模块化/循环的。

有没有一种 Python 的方式来完成所有这些映射？

---

编辑：让我们关注最困难的情况，我们有两个成对的值列表：

A        B
-----    -----
  0.0    182.5
 10.0    172.3
 20.0    161.4
 ...      ...
170.0      9.7
180.0    359.1
190.0    348.2
 ...      ...
340.0    163.6
350.0    171.8

假设列表来自两个不同的雷达，其指针未与北或其他任何方向对齐，但我们确实通过移动目标并查看每个雷达必须指向的位置来手动获取上述数据。

当雷达 A 说“我有一个 123.4 的目标！”时，我需要将雷达 B 瞄准哪里才能看到它？如果雷达 B 找到目标，我应该告诉雷达 A 指向哪里？

列表 A 环绕在最后一个元素和第一个元素之间，但列表 B 环绕在更靠近列表中间的位置。列表 A 单调增加，而列表 B 单调减少。请注意，A 上的度数的大小通常与 B 上的度数不同。

是否有一个简单的插值器可以在以下情况下正确包装：

从列表 A 插入到列表 B。

从列表 B 插值到列表 A。

可以使用两个独立的插值器实例化，一个用于每个方向。我假设线性（一阶）插值器是可以的，但我将来可能想使用高阶或样条插值。

一些测试用例：

A = 356.7, B = ?

A = 179.2, B = ?

【问题讨论】：

是否可以用一个简单的公式计算传感器读数，例如systemA = (physical*coef + offset) % 360，或者这些值是否足够非线性以使其不切实际？如果是，您可以使用代数来解决给定任何已知值的任何未知值。如果不是，那么您可能是正确的需要插值。模块化插值通常不会太糟糕，您只需要检查您之间插值的点是否相隔超过modulus/2（例如180度），表明它们之间的最短路径环绕。

你的例子没有意义。 B 的读数一直在下降，直到最后 2 个读数增加，而且它们已经小于列表顶部的读数。如果你解决了这个问题，我也许可以用我的答案展示一些示例代码。

角度数据插值的一般性注释可能会有所帮助。将数据分解为单位分量并单独对分量执行插值，然后使用扇区安全的 arctan 方法（例如 arctan2(y,x)）重新组合生成的插值，这非常有用。

【参考方案1】：

这对我有用。可能需要进行一些清理。

class InterpolatedArray(object):
    """ An array-like object that provides interpolated values between set points.
    """
    points = None
    wrap_value = None
    offset = None

    def _mod_delta(self, a, b):
        """ Perform a difference within a modular domain.
            Return a value in the range +/- wrap_value/2.
        """
        limit = self.wrap_value / 2.
        val = a - b
        if val < -limit: val += self.wrap_value
        elif val > limit: val -= self.wrap_value
        return val

    def __init__(self, points, wrap_value=None):
        """Initialization of InterpolatedArray instance.

        Parameter 'points' is a list of two-element tuples, each of which maps
        an input value to an output value.  The list does not need to be sorted.

        Optional parameter 'wrap_value' is used when the domain is closed, to
        indicate that both the input and output domains wrap.  For example, a
        table of degree values would provide a 'wrap_value' of 360.0.

        After sorting, a wrapped domain's output values must all be monotonic
        in either the positive or negative direction.

        For tables that don't wrap, attempts to interpolate values outside the
        input range cause a ValueError exception.
        """
        if wrap_value is None:
            points.sort()   # Sort in-place on first element of each tuple
        else:   # Force values to positive modular range
            points = sorted([(p[0]%wrap_value, p[1]%wrap_value) for p in points])
            # Wrapped domains must be monotonic, positive or negative
            monotonic = [points[x][1] < points[x+1][1] for x in xrange(0,len(points)-1)]
            num_pos_steps = monotonic.count(True)
            num_neg_steps = monotonic.count(False)
            if num_pos_steps > 1 and num_neg_steps > 1: # Allow 1 wrap point
                raise ValueError("Table for wrapped domains must be monotonic.")
        self.wrap_value = wrap_value
        # Pre-compute inter-value slopes
        self.x_list, self.y_list = zip(*points)
        if wrap_value is None:
            intervals = zip(self.x_list, self.x_list[1:], self.y_list, self.y_list[1:])
            self.slopes = [(y2 - y1)/(x2 - x1) for x1, x2, y1, y2 in intervals]
        else:   # Create modular slopes, including wrap element
            x_rot = list(self.x_list[1:]); x_rot.append(self.x_list[0])
            y_rot = list(self.y_list[1:]); y_rot.append(self.y_list[0])
            intervals = zip(self.x_list, x_rot, self.y_list, y_rot)
            self.slopes = [self._mod_delta(y2, y1)/self._mod_delta(x2, x1) for x1, x2, y1, y2 in intervals]

    def __getitem__(self, x):       # Works with indexing operator []
        result = None
        if self.wrap_value is None:
            if x < self.x_list[0] or x > self.x_list[-1]:
                raise ValueError('Input value out-of-range: %s'%str(x))
            i = bisect.bisect_left(self.x_list, x) - 1
            result = self.y_list[i] + self.slopes[i] * (x - self.x_list[i])
        else:
            x %= self.wrap_value
            i = bisect.bisect_left(self.x_list, x) - 1
            result = self.y_list[i] + self.slopes[i] * self._mod_delta(x, self.x_list[i])
            result %= self.wrap_value
        return result

还有一个测试：

import nose

def xfrange(start, stop, step=1.):
    """ Floating point equivalent to xrange()."""
    while start < stop:
        yield start
        start += step

# Test simple inverted mapping for non-wrapped domain
pts = [(x,-x) for x in xfrange(1.,16., 1.)]
a = InterpolatedArray(pts)
for i in xfrange(1., 15., 0.1):
    nose.tools.assert_almost_equal(a[i], -i)
# Cause expected over/under range errors
result = False  # Assume failure
try: x = a[0.5]
except ValueError: result = True
assert result
result = False
try: x = a[15.5]
except ValueError: result = True
assert result

# Test simple wrapped domain
wrap = 360.
offset = 1.234
pts = [(x,((wrap/2.) - x)) for x in xfrange(offset, wrap+offset, 10.)]
a = InterpolatedArray(pts, wrap)
for i in xfrange(0.5, wrap, 0.1):
    nose.tools.assert_almost_equal(a[i], (((wrap/2.) - i)%wrap))

【参考方案2】：

但是你可以使用线性插值。如果您的样本 A 值为例如7.75，类似于 2.5 度。如果样本 B 值为 180.35，它也类似于 2.5 度。棘手的部分是当值溢出时，如果有可能的话。只需设置一堆单元测试来检查你的算法是否有效，你应该很快就会去。

【讨论】：

所有值都在 [0,360) 范围内，其中包括 0 但不包括 360，这意味着值 360 本身从未见过，但可能会出现 359.99999999。对于一个稍微不同的问题，我必须模拟一个我们想要添加到每个系统的强度传感器，所以我有一个 0-360 范围内的广泛值，但我没有完全覆盖“背面” " 在 +/-180 度。该问题的插值器根据数据创建了一个斜率表，并使用bisect.bisect_left() 进行插值。

【参考方案3】：

第 1 部分答案：包含校准值 + 漂移值的转换表。

基本上，如果 DialA 在物理上为 0 时报告 5.7，在 5 时报告 9.7，那么我会将漂移值设置为每个读数位置之间距离的 +/- 0.25，以考虑机械和读数漂移。

第 2 部分的答案：在两个表盘上保持相同的值，同时显示预期的位置。

如果您不依赖于方向，则只需旋转输出转盘，直到它根据校准表处于正确位置。

如果您依赖于方向，则需要跟踪最后 1-2 个值以确定方向。一旦您确定了方向，您就可以沿您需要的方向移动从属转盘，直到到达目标位置。

您的校准表也应包括方向（例如正或负）。

考虑到以上两部分，您将能够补偿旋转偏移和方向翻转，并产生准确的位置和方向读数。

这里有一些代码给出了一个校准表，将产生位置和方向，这将解决显示问题并使从属表盘与主表盘匹配：

#!/usr/bin/env python

# Calibration table
# calibrations[ device ][physical position]=recorded position
calibrations=

calibrationsdrift=1.025

calibrations["WheelA"]=

calibrations["WheelA"]=
calibrations["WheelA"]["000.0"]=5.7
calibrations["WheelA"]["005.0"]=9.8
calibrations["WheelA"]["010.0"]=13.9
calibrations["WheelA"]["015.0"]=18.0

calibrations["WheelB"]=
calibrations["WheelB"]["000.0"]=182.3
calibrations["WheelB"]["005.0"]=178.4
calibrations["WheelB"]["010.0"]=174.4
calibrations["WheelB"]["015.0"]=170.3


def physicalPosition( readout , device ):
        calibration=calibrations[device]
        for physicalpos,calibratedpos in calibration.items():
                if readout < ( calibratedpos + calibrationsdrift ):
                        if readout > ( calibratedpos - calibrationsdrift ):
                                return physicalpos
        return -0


print physicalPosition( 5.8 , "WheelA")
print physicalPosition( 9.8 , "WheelA")
print physicalPosition( 10.8 , "WheelA")


def physicalDirection( lastreadout, currentreadout, device ):
        lastposition=physicalPosition( lastreadout, device)
        currentposition=physicalPosition( currentreadout, device)
        # Assumes 360 = 0, so 355 is the last position before 0
        if lastposition < currentposition:
                if lastposition == 000.0:
                        if currentposition == 355:
                                return -1
                return 1
        else:
                return -1


print physicalDirection( 5.8, 10.8, "WheelA")
print physicalDirection( 10.8, 2.8, "WheelA")

print physicalDirection( 182, 174, "WheelB")
print physicalDirection( 170, 180, "WheelB")

运行程序显示方向已正确确定，即使对于向后安装在面板/设备/等上的 WheelB 也是如此：

$ ./dials 
000.0
005.0
005.0
1
-1
1
-1

请注意，提供给函数的某些“读数”值已关闭。这是由漂移值补偿的。是否需要取决于您所连接的设备。

【讨论】：

我的参考数据每5度取一次，使用SystemA作为参考（只是一个方便的选择）。虽然物理传感器（粘贴在公共底座上的纸度计）是完全线性的，但 SystemA 和 SystemB 上的传感器都表现出超过 +/- 4 度的可变非线性，因此线性插值是绝对必要的。我正在构建一个实现来尝试将这种技术与线性插值相结合，并在我开始工作时更新我的​​问题。

【参考方案4】：

最简单的解决方案是让所有表格元素增加（或减少，视情况而定），对单个元素增加或减少 360 度以使其如此。将表格背靠背加倍，使其涵盖 0 到 360 的整个范围，即使在所有加法和减法之后也是如此。这使得简单的线性插值成为可能。然后，您可以在计算后取模 360 以将其恢复到范围内。

【讨论】：

关键的困难是知道何时/何地进行此类更正。是否有一种通用算法，即使在数据包装点也能工作？

@BobC，只需浏览表格 - 如果值从一个非常大的值跳到一个非常小的值，则从该点开始添加 360。如果它从一个小值跳到一个大值，从那个点开始减去 360。