保持 2D 插值在形状内

Posted 2023-03-27

【中文标题】保持 2D 插值在形状内

【英文标题】：Keeping 2D interpolation within shape

【发布时间】：2022-01-16 19:23:40

【问题描述】：

我有一组测量点，我想在它们之间进行插值，为此我使用 SciPy 的 griddata()：

import matplotlib.pyplot as plt 
import numpy as np
import pandas as pd 
from scipy.interpolate import griddata 

data = pd.DataFrame(
    'time': [0, 1, 2, 3, 4, 0.1, 0.9, 2, 3.05, 4, 0, 1, 2.2, 3, 3.95], 
    'force': [1, 2, 4, 9, 16, 0, 0, 0, 0, 0, -1, -2, -4, -9, -16]
)

Times, Forces = np.meshgrid(
    np.linspace(0, 4, 100),
    np.linspace(-16, 16, 100)
)

data['work'] = data['time'] * data['force']

interpolation = griddata(
    (data['time'], data['force']), 
    data['work'], 
    (Times, Forces), 
    method= 'linear'
)

fig, ax = plt.subplots() 

contour = ax.contourf(
    Times, Forces, interpolation
)
ax.scatter(data['time'], data['force'])

fig.show()

我的问题是，我的测量点已经遵循物理可能性的边界，但是插值仍然会针对每个可跨区域进行插值，包括那些无法测量的区域。

如何将插值或至少插值图限制在外部点的形状“范围内”？不幸的是，“时间”测量值有小的偏差。

任何建议将不胜感激！

【问题讨论】：

你试过ax.tricontourf(data['time'], data['force'], data['work'],...)吗？这将使用给定坐标的凸包的三角剖分。

@JohanC 我没有，但恐怕它也有同样的问题。但无论如何我可能会使用它，因为它是我迄今为止所写的更优雅的版本，所以绝对感谢你！ :) 但也许我可以创建一个可以覆盖出界区域的面具......

它可能会在work 值np.nan 的边界区域中添加点。

【参考方案1】：

首先，您需要获取所有点的凹多边形。其次，使用多边形来裁剪轮廓填充。虽然这可能有点复杂，但一些有用的包可以帮助完成这些任务。

下面是代码。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.interpolate import griddata
import matplotlib.path as mpath

import alphashape
from descartes import PolygonPatch

Path = mpath.Path

data = pd.DataFrame(
    'time': [0, 1, 2, 3, 4, 0.1, 0.9, 2, 3.05, 4, 0, 1, 2.2, 3, 3.95],
    'force': [1, 2, 4, 9, 16, 0, 0, 0, 0, 0, -1, -2, -4, -9, -16]
)

Times, Forces = np.meshgrid(
    np.linspace(0, 4, 100),
    np.linspace(-16, 16, 100)
)

data['work'] = data['time'] * data['force']

interpolation = griddata(
    (data['time'], data['force']),
    data['work'],
    (Times, Forces),
    method= 'linear')

fig, ax = plt.subplots()
plt.xlim(-.3,4.3)
plt.ylim(-18,18)

#contour = ax.contourf( Times, Forces, interpolation)
contour = ax.tricontourf( data['time'], data['force'], data['work'])
ax.scatter(data['time'], data['force'])


x = data['time']
y =  data['force']

points = np.vstack([x, y]).T
alpha = 0.95 * alphashape.optimizealpha(points)
hull = alphashape.alphashape(points, alpha)
hull_pts = hull.exterior.coords.xy

ax.scatter(hull_pts[0], hull_pts[1], color='red')
ax.add_patch(PolygonPatch(hull, fill=False, color='red'))

plt.savefig("clip_before.png")
#make clip path
vertices = []
codes = []
xpts,ypts = hull_pts

## convert polygon to path for cliping contour fill
for ix,iy in zip(xpts,ypts):
    vertices.append((ix,iy))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (len(xpts) -2)
codes += [Path.CLOSEPOLY]
clip = Path(vertices, codes)

for collection in contour.collections:
    collection.set_clip_path(clip,transform=ax.transData)
plt.savefig("clip_after.png")
plt.show()

这是输出数字。

剪裁前。

剪裁后。

【讨论】：

真的很优雅，非常感谢！但是有一个问题：创建的多边形的形状受 alpha 的影响，对吗？为什么是 0.95，多边形也可以定义为通过更多点吗？