平滑二维图形

Question

我有一些含糊不清的矩形二维图形需要平滑处理。以下是一个简化的例子：

fig, ax1 = plt.subplots(1,1, figsize=(3,3))
xs1 = [-0.25,  -0.625, -0.125, -1.25, -1.125, -1.25, 0.875, 1.0, 1.0, 0.5, 1.0, 0.625, -0.25]
ys1 = [1.25, 1.375, 1.5, 1.625, 1.75, 1.875, 1.875, 1.75, 1.625, 1.5, 1.375, 1.25, 1.25]

ax1.plot(xs1, ys1)
ax1.set_ylim(0.5,2.5)
ax1.set_xlim(-2,2) ;

我尝试了scipy.interpolate.RectBivariateSpline，但这显然需要所有点的数据（例如热图），而scipy.interpolate.interp1d则合理地要求生成1D平滑版本。

有什么适当的方法来平滑处理这个问题？

编辑以更好地修订/解释我的目标。 我不需要线穿过所有点; 实际上，我希望它们不要穿过所有点，因为有明显的异常值点应该与邻居平均或采用类似的方法。 我在上面包含了一个粗略的手动草图，作为我所想象的开始。

- iayork

Answer 1

Chaikin's corner cutting algorithm是您可能的理想方法。对于一个由顶点P0、P1…P(N-1)定义的多边形，角切割算法将为由P(i)和P(i+1)定义的每个线段生成2个新的顶点:

Q(i) = (3/4)P(i) + (1/4)P(i+1)
R(i) = (1/4)P(i) + (3/4)P(i+1)

因此，您的新多边形将有2N个顶点。然后，您可以再次在新的多边形上应用角切割，并反复进行，直到达到期望的分辨率。结果将是一个具有许多顶点的多边形，但它们在显示时看起来很平滑。可以证明，使用这种角切割方法产生的曲线会收敛成二次B样条曲线。该方法的优点是所得到的曲线永远不会过度调整。下面的图片将为您提供更好的关于此算法的概念（图片来自此链接）：

原始多边形

应用一次角切割

再应用一次角切割

有关Chaikin的角切割算法的更多详细信息，请参见此链接。

- fang

---

这是该算法的Python 2.7.X实现：https://dev59.com/kKbja4cB1Zd3GeqPgG-w - Comrade Che

Answer 2

实际上，您可以使用scipy.interpolate.inter1d来实现此目的。您需要将多边形的x和y组件分别视为不同的系列。

以下是一个快速示例，使用正方形：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

x = [0, 1, 1, 0, 0]
y = [0, 0, 1, 1, 0]

t = np.arange(len(x))
ti = np.linspace(0, t.max(), 10 * t.size)

xi = interp1d(t, x, kind='cubic')(ti)
yi = interp1d(t, y, kind='cubic')(ti)

fig, ax = plt.subplots()
ax.plot(xi, yi)
ax.plot(x, y)
ax.margins(0.05)
plt.show()

然而，如您所见，这在0,0处会出现一些问题。

这是由于样条段不仅依赖于两个点。我们所插值的第一个和最后一个点并没有以我们期望的方式“连接”。我们可以通过在x和y序列中添加倒数第二个和第二个点来解决这个问题，并为样条线端点提供“包裹式”边界条件。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

x = [0, 1, 1, 0, 0]
y = [0, 0, 1, 1, 0]

# Pad the x and y series so it "wraps around".
# Note that if x and y are numpy arrays, you'll need to
# use np.r_ or np.concatenate instead of addition!
orig_len = len(x)
x = x[-3:-1] + x + x[1:3]
y = y[-3:-1] + y + y[1:3]

t = np.arange(len(x))
ti = np.linspace(2, orig_len + 1, 10 * orig_len)

xi = interp1d(t, x, kind='cubic')(ti)
yi = interp1d(t, y, kind='cubic')(ti)

fig, ax = plt.subplots()
ax.plot(xi, yi)
ax.plot(x, y)
ax.margins(0.05)
plt.show()

只是为了展示使用您的数据时的样子：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

x = [-0.25, -0.625, -0.125, -1.25, -1.125, -1.25,
     0.875, 1.0, 1.0, 0.5, 1.0, 0.625, -0.25]
y = [1.25, 1.375, 1.5, 1.625, 1.75, 1.875, 1.875,
     1.75, 1.625, 1.5, 1.375, 1.25, 1.25]

# Pad the x and y series so it "wraps around".
# Note that if x and y are numpy arrays, you'll need to
# use np.r_ or np.concatenate instead of addition!
orig_len = len(x)
x = x[-3:-1] + x + x[1:3]
y = y[-3:-1] + y + y[1:3]

t = np.arange(len(x))
ti = np.linspace(2, orig_len + 1, 10 * orig_len)

xi = interp1d(t, x, kind='cubic')(ti)
yi = interp1d(t, y, kind='cubic')(ti)

fig, ax = plt.subplots()
ax.plot(xi, yi)
ax.plot(x, y)
ax.margins(0.05)
plt.show()

请注意，这种方法会导致相当多的“超调”。这是由于三次样条插值造成的。@pythonstarter的建议也是处理它的另一种好方法，但贝塞尔曲线将遇到相同的问题（它们在数学上基本上是等价的，只是控制点的定义方式不同）。还有许多其他处理平滑的方法，包括专门用于平滑多边形的方法（例如指数核多项式逼近（PAEK））。我从未尝试过实现PAEK，所以我不确定它的复杂程度。如果您需要更强大的平滑功能，可以尝试查找PAEK或其他类似的方法。

- Joe Kington

---

这个方案是可行的，并且符合我最初解释的描述。然而，事实证明我真正需要的略有不同，@fang提供的Chaikin拐角削减算法更接近我想要的。 - iayork

这解决了我的贝塞尔曲线问题！我卡了整整一个晚上，试图弄清楚如何使循环边界条件起作用，但一无所获... - Flying_Banana

Answer 3

这更像是一条评论而不是答案，但也许你可以尝试将该多边形定义为贝塞尔曲线。代码相当简单，我相信您熟悉这些曲线的工作原理。在这种情况下，该曲线将成为控制多边形。但并非所有都完美：首先，它实际上不是该多边形的“平滑版本”，而是一条曲线；另一件事是，曲线的阶数越高，它看起来就越不像控制多边形。我的意思是，也许您应该尝试使用数学工具来解决这个问题，而不是用编程技巧来平滑多边形。

- python starter