检测阵列的抛物面形状部分

Question

我试图在^{}中检测抛物面形状的点。你知道吗

我有一个数据集，使用维纳过滤器（使用包NIFTy4）对其进行平滑处理：它被保存为.csv，可以在here on Google Drive或here on PasteBin中找到。你知道吗

我要寻找的是一种方法来识别数组中描述向上打开的抛物面的部分。你知道吗

在提供的样本中，我想检测5个这样的形状。重点主要是形状的宽度（开始和结束），而不是实际的最小位置。你知道吗

提前谢谢

Answer 1

据我所知，曲线是平滑的，无噪声的，因此抛物线顶点只对应于局部极小值。你知道吗

您没有指定如何定义开始和结束，因此我将假设“直到拐点”。要定位这些，可以计算一阶差来估计坡度，并检测两侧绝对坡度的第一个最大值。你知道吗

不要期望在这些极限上有很好的精度，这是一个不适定的问题。你知道吗

Answer 2

使用numpy的解决方案

对于这个解决方案，我将使用numpy：

import numpy as np

数据集

创建数据集

OPPhteven很好地提供了要处理的数据集，但是由于到数据集的链接趋于消失，我还创建了一个函数来生成类似的曲线。你知道吗

def polyval_points(n):
    """Return random points generated from polyval."""
    x = np.linspace(-20, 20, num=n)
    coef = np.random.normal(-3, 3, size=(5))
    y = np.polyval(coef, x) * np.sin(x)
    return x, y

加载数据集

def dataset_points():
    """Return points loaded from dataset."""
    y = np.loadtxt("data.csv", delimiter=',')
    x = np.linspace(0, 1, num=len(y))
    return x, y

阐述要点

斜率卷积

因为这些点是离散的，所以我们必须平均斜率。其中一种方法是通过统一内核。你知道吗

def convolute_slopes(y, k=3):
    """Return slopes convoluted with an uniform kernel of size k."""
    d2y = np.gradient(np.gradient(y))
    return np.convolve(d2y, np.ones(k)/k, mode="same")

获得抛物面

现在，我们可以计算卷积斜率，确定它在何处切换方向，并选择那些平均斜率大于绝对平均值乘以一个描述抛物面“斜率”的系数的间隔。你知道吗

def get_paraboloids(x, y, c=0.2):
    """Return list of points (x,y) that are part of paraboloids in given set.
        x: np.ndarray of floats
        y: np.ndarray of floats
        c: slopyness coefficient
    """
    slopes = convolute_slopes(y)
    mean = np.mean(np.abs(slopes))
    w = np.where(np.diff(slopes > 0) > 0)[0] + 1
    w = np.insert(w, [0, len(w)], [0, len(x)])
    return [(x[lower:upper], y[lower:upper])
            for lower, upper in zip(w[:-1], w[1:])
            if np.mean(slopes[lower:upper]) > mean * c]

如何使用这个和可视化

首先，我们加载数据集并生成更多数据：

datasets = [dataset_points(), polyval_points(10000), polyval_points(10000)]

然后，迭代每个数据集：

from matplotlib import pyplot as plt

plt.figure(figsize=(10 * len(datasets), 10))
for i, points in enumerate(datasets):
    x, y = points
    plt.subplot(1, len(datasets), i + 1)
    plt.plot(x, y, linewidth=1)
    for gx, gy in get_paraboloids(x, y):
        plt.plot(gx, gy, linewidth=3)
plt.show()

结果