如何估计有噪口罩的分段平滑度？

Question

我有一个源自分段平滑曲线的二进制掩码。每个基础曲线线段都是平滑的，但是线段之间的连接可能是平滑的，也可能不是平滑的。遮罩很吵，因此它可能在基础曲线周围包含多个像素。下图显示了此类输入的示例：在没有任何其他先验知识的情况下，鉴于此输入，我想估计与基础曲线的拟合，从而能够提供平滑和非平滑的连接。

我正在python中工作，所以我尝试了几种可用的方法，例如numpy的多项式拟合，scipy的样条平滑和pyqt-fit的非参数回归，但是无法获得所需的输出。这是一个代码示例：

    from imageio import imread
    import random
    import numpy as np
    from scipy.interpolate import UnivariateSpline
    import pyqt_fit.nonparam_regression as smooth
    from pyqt_fit import npr_methods
    from matplotlib import pyplot as plt

    mask = imread(r'C:\mask.bmp')
    ys, xs = np.where(mask)
    # add tiny noise and sort - silly way to comply to UnivariateSpline's requirement of "x must be strictly increasing"
    xs = np.array([x + random.random() * 1e-4 for x in xs])
    sorter = np.argsort(xs)
    xs = xs[sorter]
    ys = ys[sorter]
    # polynomial fit
    p = np.poly1d(np.polyfit(xs, ys, 5))
    # spline smoothing
    spl = UnivariateSpline(xs, ys, k=3, s=1e9)
    # non-parameteric regression
    k = smooth.NonParamRegression(xs, ys, method=npr_methods.LocalPolynomialKernel(q=3))
    k.fit()

    plt.figure()
    plt.imshow(mask, cmap='gray')
    linexs = np.array(range(mask.shape[1]))
    plt.plot(linexs, k(linexs), 'y', lw=1)
    plt.plot(linexs, spl(linexs), 'g', lw=1)
    plt.plot(linexs, p(linexs), 'b', lw=1)
    plt.show()

对于此示例中显示的参数，这些配合均无法捕获左侧的非平滑连接，也无法为右侧的“尾部”提供良好的配合：预期结果的行为应类似于下图中的红色曲线，在这里我希望位置1处的曲线不平滑，位置2处的曲线平滑。

我很乐意引用合适的算法。如果还有python实现，那将是一个加分。

Answer 1

这可能接近您想要的。诀窍是首先skeletonize蒙版（请参见https://scikit-image.org/docs/dev/auto_examples/edges/plot_skeleton.html）以提取要估计曲线的点。

from skimage.morphology import skeletonize
import matplotlib.pyplot as plt
import cv2
import numpy as np
from scipy.interpolate import interp1d

img = cv2.imread("./data/bAiew.png", 0)
img = cv2.medianBlur(img, 11)
img = img//255

# Sparse skeleton
skeleton = skeletonize(img, method='lee')

# Into dense
contours, _ = cv2.findContours(skeleton, 0,cv2.CHAIN_APPROX_NONE)
arch = contours[0]
x,y = arch[...,0].squeeze(), arch[...,1].squeeze()

# Fitting a curve
xx, yy = x[0::15], y[0::15] #<- sample every 15th element to see that the interpolate really works
f = interp1d(xx, yy)

plt.figure(figsize=(16,9))
plt.subplot(221)
plt.imshow(img)
plt.title('original mask')

plt.subplot(222)
plt.imshow(skeleton)
plt.title('skeleton')

plt.subplot(223)
plt.scatter(xx,yy)
plt.ylim([img.shape[0],0])
plt.xlim([0, img.shape[1]])
plt.title('sampled points')

plt.subplot(224)
plt.plot(xx,f(xx))
plt.ylim([img.shape[0],0])
plt.xlim([0, img.shape[1]])
plt.title('interpolate')
plt.show()

输出：