python - 检测阵列的抛物面形部分

Question

我正在尝试检测numpy.ndarray.

我有一个数据集，使用 Wiener 过滤器（使用包 NIFTy4）对其进行平滑处理：它被保存为.csv并且可以在 Google Drive或PasteBin 上找到。

我正在寻找的是一种方法来识别我的阵列中描述向上打开的抛物面的部分。

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

提前致谢

标签： pythonnumpypattern-matching

Answer 1

解决方案使用numpy

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

import numpy as np

数据集

创建数据集

OP Phteven非常友好地提供了可以使用的数据集，但由于与数据集的链接往往会消失，我还创建了一个生成类似曲线的函数。

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()

结果