【问题标题】:Interpolation and Extrapolation of Randomly Scattered data to Uniform Grid in 3D随机分散数据的插值和外推到 3D 均匀网格
【发布时间】:2014-03-26 11:07:28
【问题描述】:

我有一个 256 x 256 x 32 的规则间隔点网格,分布在 x、y 和 z 上,并带有关联的变量“a”。我在更有限的 x、y、z 空间中也有一组随机分散的点,并带有关联的变量“b”。我本质上想要做的是将我的随机数据内插和外推到与“a”立方体匹配的规则间隔网格中,如下所示:

到目前为止,我已经使用 scipy 的 griddata 来实现插值,这似乎工作正常,但它无法处理外插(据我所知)并且输出急剧截断为“nan”值。在研究这个问题时,我遇到了几个人第二次使用 griddata 使用“最近”作为插值方法来填充“nan”值。我试过这个,但结果似乎不可靠。如果我使用具有“线性”模式的 fill_Value,可以获得更合适的外观结果,但目前它更像是一个软糖,因为 fill_Value 必须是一个常数。

我注意到 MATLAB 有一个 ScatteredInterpolant 类,它似乎可以满足我的要求,但我无法在 Python 中找到等效的类,也无法弄清楚如何在 3D 中有效地实现这样的例程。非常感谢任何帮助。

我用于插值的代码如下:

x, y, z, b = np.loadtxt(scatteredfile, unpack = True)

# Create cube to match aCube dimensions
xi = np.linspace(-xmax_aCube, xmax_aCube, 256)
yi = np.linspace(-ymax_aCube, ymax_aCube, 256)
zi = np.linspace(zmin_aCube, zmax_aCube, 32)

# Interpolate scattered points
X, Y, Z = np.meshgrid(xi, yi, zi)
bCube = griddata((x, y, z), b, (X, Y, Z), method = 'linear')    

【问题讨论】:

  • 是空间计数数据中的点还是它们与相应的测量坐标f(x,y,z)?你确定这是插值和外推吗?看起来您可能希望将这些点放入您构建的 3-d 箱中。
  • 如果您的 b 数据表现良好,您可以向分散的云添加一些均匀分布的合成数据点以强制执行您想要的行为 - 但这完全取决于您的用例(可能有理由简单地显示数据)...

标签: python scipy interpolation extrapolation


【解决方案1】:

这个讨论适用于任何维度。对于您的 3D 案例,让我们先谈谈计算几何,以了解为什么部分区域从 griddata 给出 NaN

你体积中的散点组成一个凸包;具有以下属性的几何形状:

  • 表面总是凸的(顾名思义)
  • 在不违反凸度的情况下,形状的体积尽可能小
  • 表面(在 3d 中)被三角剖分并闭合

不太正式,凸包 (which you can compute easily with scipy) 就像在框架上拉伸气球,其中框架角是分散集群的最外点。

在常规网格位置内部,您被已知点包围。您可以对这些位置进行插值。在它之外,您必须进行推断。

外推很难。没有关于如何做的一般规则......它是针对特定问题的。在那个区域,像griddata 这样的算法选择 返回NaN - 这是告知科学家他/她必须选择合理的推断方式的最安全方式。

让我们通过一些方法来做到这一点。

1。 [最糟糕] 搞砸了

在外壳外部分配一些标量值。 In the numpy docs 你会看到这是通过以下方式完成的: s = 平均值(b) bCube = griddata((x, y, z), b, (X, Y, Z), method = 'linear', fill_value=s)

缺点:这会在船体边界处的插值场中产生明显的不连续性,严重偏离平均标量场值并且不尊重数据的函数形式。

2。 [NEXT WORST] “混合搞砸了”

假设您在域的角落应用了一些值。这可能是与散点相关的标量场的平均值。

对不起,这是伪代码,因为我根本不使用 numpy,但它可能会很清楚

# With a unit cube, and selected scalar value
x, y, z, b = np.loadtxt(scatteredfile, unpack = True)
s = mean(b)
x.append([0 0 0 0 1 1 1 1])
y.append([0 0 1 1 0 0 1 1])
z.append([0 1 0 1 0 1 0 1])
b.append([s s s s s s s s])
# drop in the rest of your code

缺点:这会在船体边界处的插值场的梯度中产生明显的不连续性,相当大地偏离平均标量场值并且不尊重函数数据的形式。

3。 [STILL PRETTY BAD] 最近的邻居

对于每个常规 NaN 点,找到最近的非 NaN 并分配该值。这是有效且稳定的,但很粗糙,因为您的字段最终可能会出现图案化特征(例如从船体向外辐射的条纹或光束),通常在视觉上不吸引人,或者更糟糕的是,在数据平滑度方面无法接受

根据数据的密度,您可以使用最近的分散数据点而不是最近的非 NaN 常规点。这可以通过(再次,伪代码)简单地完成:

bCube = griddata((x, y, z), b, (X, Y, Z), method = 'linear', fill_value=nan)
bCubeNearest = griddata((x, y, z), b, (X, Y, Z), method = 'nearest')
indicesMask = isNan(bCube)
# Use nearest interpolation outside the hull, keeping linear interpolation inside.
bCube(indicesMask) = bCubeNearest(indicesMask)

使用 MATLAB 的基于 delaunay 的方法将揭示更强大的方法来实现类似的单行,但 numpy 在这里看起来有点有限。

4。 [不总是很糟糕] 自然加权

抱歉这部分的解释不好,我从来没有写过算法,但我相信对自然邻域技术的一些研究会让你走得更远

使用带有某个参数D 的距离加权函数,该参数可能与您的盒子长度相似或两倍(比如)。你可以调整。对于每个 NaN 位置,计算到每个散点的距离。

# Don't do it this way for anything but small matrices - this is O(NM) 
# and it can be done much more effectively (e.g. MATLAB has a quick 
# natural weighting option), but for illustrative purposes:
for each NaN point 1:N
    for each scattered point 1:M
        calculate a basis function using inverse distance from NaN to point, normalised on D, and store in a [1 x M] vector of weights
Multiply weights by the b value, summate and divide by M   

您基本上希望得到一个函数,该函数在距船体距离 D 处平滑地达到 B 的平均强度,但与边界处的船体重合。远离边界,它在最近的点上的权重最大。

优点:非常稳定且相当连续。由于加权,在单个数据点上比最近邻更能抵御噪声。

5。 [HEROIC ROCKSTAR] 函数形式假设

你对物理学了解多少?假设一个函数形式代表你期望物理学做的事情,然后对该形式进行最小二乘(或某种等效)拟合到分散的数据。使用函数来稳定外推。

一些可以帮助你构建函数的好主意:

  • 您期望对称性还是周期性?
  • b 是矢量场的一个分量,它具有诸如零散度之类的特性吗?
  • 方向性:您是否希望所有角都相同?或者可能是一个方向的线性变化?
  • 字段 b 是否在某个时间点 - 也许可以使用平滑的测量时间序列来得出基本函数?
  • 是否已经存在已知的形式,例如高斯或二次?

一些例子:

  • b 表示通过体积的激光束的强度。您希望入口侧名义上与出口相同,其他四个边界为零强度。强度将具有同心高斯分布。

  • b 是不可压缩流体中速度场的一个分量。流体必须是无发散的,因此在 NaN 区域产生的任何场也必须是无发散的,因此您可以应用此条件。

  • b 代表房间的温度。您预计顶部的温度会更高,因为热空气会上升。

  • b 代表机翼上的升力,在三个独立变量上进行了测试。您可以轻松地在失速时查看电梯,因此确切地知道它在空间的某些部分会是什么。

优点/缺点:做对了,它会很棒。弄错了,尤其是非线性函数形式,它会非常出错,并可能导致非常不稳定的结果。

健康警告你不能假设一个函数形式,得到漂亮的结果,然后用它们来证明函数形式是正确的。那只是糟糕的科学。表单必须是行为良好且已知独立于您的数据分析的东西。

【讨论】:

    【解决方案2】:

    如果您的点散布点非常符合立方体形状,一种方法可能是使用 griddata 插入到适合您的点云的常规数据网格(因此避免使用 nan),然后使用此常规网格值作为interpn 的输入,这确实有助于线性外推(但需要常规网格作为输入)。

    这样,您可以像以前一样将griddata 用于点散布的凸包内的所有点,并且可以使用interpn 来估计以nans 形式返回的点。

    这远非完美,但我认为它更接近于实现您的目标。

    优点:

    • 避免明显的不连续性。
    • 无需了解函数形式即可在数据集边缘捕获基本线性趋势。
    • 尊重数据中的不对称性(例如,不倾向于远距离的总体均值,因此数据集的一侧在远距离时可能比另一侧具有更大的值。)

    缺点:

    • 此方法的有效性在很大程度上取决于您可以在初始散布点的凸包内容纳多大的立方体。如果您的数据是尖峰/不规则且不规则的,则凸包边缘上的甚至点可能已被外推到距嵌套立方体边缘很远的距离,从而导致错误,因为外推不会考虑位于较近的数据点立方体外。
    • 线性外推会受到数据中噪声的严重影响 在点云的边缘。
    • 进行两组插值的计算成本。

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 2016-06-17
      • 1970-01-01
      • 2021-07-20
      • 1970-01-01
      • 1970-01-01
      • 2011-08-08
      • 1970-01-01
      • 2020-04-04
      相关资源
      最近更新 更多